From 092944888a35636d341525d1f0f8b1a812a2b3ea Mon Sep 17 00:00:00 2001 From: Ramiro Torres-Escobedo Date: Thu, 15 Sep 2022 15:17:09 -0500 Subject: [PATCH 1/6] adding functionality for n_transits --- hawc_hal/HAL.py | 181 +++++++------------- hawc_hal/maptree/from_hdf5_file.py | 77 +++++---- hawc_hal/maptree/from_root_file.py | 266 ++--------------------------- hawc_hal/maptree/map_tree.py | 74 ++++---- hawc_hal/psf_fast/psf_wrapper.py | 44 +---- hawc_hal/response/response.py | 197 +++++++-------------- hawc_hal/response/response_bin.py | 127 +------------- hawc_hal/tests/conftest.py | 87 +++++----- hawc_hal/tests/test_n_transits.py | 33 ++++ 9 files changed, 327 insertions(+), 759 deletions(-) create mode 100644 hawc_hal/tests/test_n_transits.py diff --git a/hawc_hal/HAL.py b/hawc_hal/HAL.py index 03de58e..22ed34f 100644 --- a/hawc_hal/HAL.py +++ b/hawc_hal/HAL.py @@ -25,12 +25,17 @@ from threeML.utils.statistics.gammaln import logfactorial from tqdm.auto import tqdm -from hawc_hal.convolved_source import (ConvolvedExtendedSource2D, - ConvolvedExtendedSource3D, - ConvolvedPointSource, - ConvolvedSourcesContainer) -from hawc_hal.healpix_handling import (FlatSkyToHealpixTransform, - SparseHealpix, get_gnomonic_projection) +from hawc_hal.convolved_source import ( + ConvolvedExtendedSource2D, + ConvolvedExtendedSource3D, + ConvolvedPointSource, + ConvolvedSourcesContainer, +) +from hawc_hal.healpix_handling import ( + FlatSkyToHealpixTransform, + SparseHealpix, + get_gnomonic_projection, +) 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 @@ -53,19 +58,34 @@ 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=None, + ): # Store ROI self._roi = roi + # optionally specify n_transits + if set_transits is not None: + log.info(f"Setting transits to {set_transits}") + n_transits = set_transits + + else: + n_transits = None + log.info("Using transits contained in maptree") + # 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 - ) + self._flat_sky_projection = self._roi.get_flat_sky_projection(self.flat_sky_pixels_size) # Read map tree (data) - self._maptree = map_tree_factory(maptree, roi=self._roi) + self._maptree = map_tree_factory(maptree, roi=self._roi, n_transits=n_transits) # Read detector response_file self._response = hawc_response_factory(response_file) @@ -186,9 +206,7 @@ def psf_integration_method(self, mode): def _setup_psf_convolutors(self): - central_response_bins = self._response.get_response_dec_bin( - self._roi.ra_dec_center[1] - ) + central_response_bins = self._response.get_response_dec_bin(self._roi.ra_dec_center[1]) self._psf_convolutors = collections.OrderedDict() for bin_id in central_response_bins: @@ -261,9 +279,7 @@ def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=Non for this_bin in range(bin_id_min, bin_id_max + 1): this_bin = str(this_bin) if this_bin not in self._all_planes: - raise ValueError( - f"Bin {this_bin} is not contained in this maptree." - ) + raise ValueError(f"Bin {this_bin} is not contained in this maptree.") self._active_planes.append(this_bin) @@ -281,9 +297,7 @@ def set_active_measurements(self, bin_id_min=None, bin_id_max=None, bin_list=Non # if not this_bin in self._all_planes: if this_bin not in self._all_planes: - raise ValueError( - f"Bin {this_bin} is not contained in this maptree." - ) + raise ValueError(f"Bin {this_bin} is not contained in this maptree.") self._active_planes.append(this_bin) @@ -438,9 +452,7 @@ def get_excess_background(self, ra: float, dec: float, radius: float): # each radial bin data = data_analysis_bin.observation_map.as_partial() bkg = data_analysis_bin.background_map.as_partial() - mdl = self._get_model_map( - energy_id, n_point_sources, n_ext_sources - ).as_partial() + mdl = self._get_model_map(energy_id, n_point_sources, n_ext_sources).as_partial() # select counts only from the pixels within specifid distance from # origin of radial profile @@ -507,10 +519,7 @@ def get_radial_profile( # The area of each ring is then given by the difference between two # subsequent circe areas. area = np.array( - [ - self.get_excess_background(ra, dec, r + offset * delta_r)[0] - for r in radii - ] + [self.get_excess_background(ra, dec, r + offset * delta_r)[0] for r in radii] ) temp = area[1:] - area[:-1] @@ -518,10 +527,7 @@ def get_radial_profile( # signals signal = np.array( - [ - self.get_excess_background(ra, dec, r + offset * delta_r)[1] - for r in radii - ] + [self.get_excess_background(ra, dec, r + offset * delta_r)[1] for r in radii] ) temp = signal[1:] - signal[:-1] @@ -529,10 +535,7 @@ def get_radial_profile( # backgrounds bkg = np.array( - [ - self.get_excess_background(ra, dec, r + offset * delta_r)[2] - for r in radii - ] + [self.get_excess_background(ra, dec, r + offset * delta_r)[2] for r in radii] ) temp = bkg[1:] - bkg[:-1] @@ -543,10 +546,7 @@ def get_radial_profile( # model # convert 'top hat' excess into 'ring' excesses. model = np.array( - [ - self.get_excess_background(ra, dec, r + offset * delta_r)[3] - for r in radii - ] + [self.get_excess_background(ra, dec, r + offset * delta_r)[3] for r in radii] ) temp = model[1:] - model[:-1] @@ -557,10 +557,7 @@ def get_radial_profile( self.set_model(model_to_subtract) model_subtract = np.array( - [ - self.get_excess_background(ra, dec, r + offset * delta_r)[3] - for r in radii - ] + [self.get_excess_background(ra, dec, r + offset * delta_r)[3] for r in radii] ) temp = model_subtract[1:] - model_subtract[:-1] @@ -580,17 +577,11 @@ def get_radial_profile( # them to the data later. Weight is normalized (sum of weights over # the bins = 1). - total_excess = np.array(self.get_excess_background(ra, dec, max_radius)[1])[ - good_planes - ] + total_excess = np.array(self.get_excess_background(ra, dec, max_radius)[1])[good_planes] - total_bkg = np.array(self.get_excess_background(ra, dec, max_radius)[2])[ - good_planes - ] + total_bkg = np.array(self.get_excess_background(ra, dec, max_radius)[2])[good_planes] - total_model = np.array(self.get_excess_background(ra, dec, max_radius)[3])[ - good_planes - ] + total_model = np.array(self.get_excess_background(ra, dec, max_radius)[3])[good_planes] w = np.divide(total_model, total_bkg) weight = np.array([w / np.sum(w) for _ in radii]) @@ -644,13 +635,7 @@ def plot_radial_profile( values for easy retrieval """ - ( - radii, - excess_model, - excess_data, - excess_error, - plane_ids, - ) = self.get_radial_profile( + (radii, excess_model, excess_data, excess_error, plane_ids,) = self.get_radial_profile( ra, dec, active_planes, @@ -682,9 +667,7 @@ def plot_radial_profile( plt.plot(radii, excess_model, color="red", label="Model") - plt.legend( - bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16 - ) + plt.legend(bbox_to_anchor=(1.0, 1.0), loc="upper right", numpoints=1, fontsize=16) plt.axhline(0, color="deepskyblue", linestyle="--") x_limits = [0, max_radius] @@ -794,9 +777,7 @@ def display_spectrum(self): yerr = [yerr_high, yerr_low] - return self._plot_spectrum( - net_counts, yerr, model_only, residuals, residuals_err - ) + return self._plot_spectrum(net_counts, yerr, model_only, residuals, residuals_err) def _plot_spectrum(self, net_counts, yerr, model_only, residuals, residuals_err): @@ -870,9 +851,7 @@ def get_log_like(self, individual_bins=False, return_null=False): bkg_renorm = list(self._nuisance_parameters.values())[0].value obs = data_analysis_bin.observation_map.as_partial() # type: np.array - bkg = ( - data_analysis_bin.background_map.as_partial() * bkg_renorm - ) # type: np.array + bkg = data_analysis_bin.background_map.as_partial() * bkg_renorm # type: np.array this_pseudo_log_like = log_likelihood(obs, bkg, this_model_map_hpx) @@ -965,9 +944,7 @@ def get_simulated_dataset(self, name): # Active plane. Generate new data expectation = self._clone[1][bin_id] - new_data = np.random.poisson( - expectation, size=(1, expectation.shape[0]) - ).flatten() + new_data = np.random.poisson(expectation, size=(1, expectation.shape[0])).flatten() # Substitute data data_analysis_bin.observation_map.set_new_values(new_data) @@ -977,18 +954,16 @@ def get_simulated_dataset(self, name): # Adjust the name of the nuisance parameter old_name = list(self._clone[0]._nuisance_parameters.keys())[0] new_name = old_name.replace(self.name, name) - self._clone[0]._nuisance_parameters[new_name] = self._clone[ - 0 - ]._nuisance_parameters.pop(old_name) + self._clone[0]._nuisance_parameters[new_name] = self._clone[0]._nuisance_parameters.pop( + old_name + ) # Recompute biases self._clone[0]._compute_likelihood_biases() return self._clone[0] - def _get_expectation( - self, data_analysis_bin, energy_bin_id, n_point_sources, n_ext_sources - ): + def _get_expectation(self, data_analysis_bin, energy_bin_id, n_point_sources, n_ext_sources): # Compute the expectation from the model @@ -1004,9 +979,7 @@ def _get_expectation( psf_integration_method=self._psf_integration_method, ) - expectation_from_this_source = ( - expectation_per_transit * data_analysis_bin.n_transits - ) + expectation_from_this_source = expectation_per_transit * data_analysis_bin.n_transits if this_model_map is None: @@ -1045,18 +1018,14 @@ def _get_expectation( # Only extended sources this_model_map = ( - self._psf_convolutors[energy_bin_id].extended_source_image( - this_ext_model_map - ) + self._psf_convolutors[energy_bin_id].extended_source_image(this_ext_model_map) * data_analysis_bin.n_transits ) else: this_model_map += ( - self._psf_convolutors[energy_bin_id].extended_source_image( - this_ext_model_map - ) + self._psf_convolutors[energy_bin_id].extended_source_image(this_ext_model_map) * data_analysis_bin.n_transits ) @@ -1066,18 +1035,14 @@ def _get_expectation( # First divide for the pixel area because we need to interpolate brightness # this_model_map = old_div(this_model_map, self._flat_sky_projection.project_plane_pixel_area) - this_model_map = ( - this_model_map / self._flat_sky_projection.project_plane_pixel_area - ) + this_model_map = this_model_map / self._flat_sky_projection.project_plane_pixel_area this_model_map_hpx = self._flat_sky_to_healpix_transform[energy_bin_id]( this_model_map, fill_value=0.0 ) # Now multiply by the pixel area of the new map to go back to flux - this_model_map_hpx *= hp.nside2pixarea( - data_analysis_bin.nside, degrees=True - ) + this_model_map_hpx *= hp.nside2pixarea(data_analysis_bin.nside, degrees=True) else: @@ -1152,9 +1117,7 @@ def display_fit(self, smoothing_kernel_sigma=0.1, display_colorbar=False): this_ra, this_dec = self._roi.ra_dec_center # Make a full healpix map for a second - whole_map = self._get_model_map( - plane_id, n_point_sources, n_ext_sources - ).as_dense() + whole_map = self._get_model_map(plane_id, n_point_sources, n_ext_sources).as_dense() # Healpix uses longitude between -180 and 180, while R.A. is between 0 and 360. We need to fix that: longitude = ra_to_longitude(this_ra) @@ -1163,9 +1126,7 @@ def display_fit(self, smoothing_kernel_sigma=0.1, display_colorbar=False): latitude = this_dec # Background and excess maps - bkg_subtracted, _, background_map = self._get_excess( - data_analysis_bin, all_maps=True - ) + bkg_subtracted, _, background_map = self._get_excess(data_analysis_bin, all_maps=True) # Make all the projections: model, excess, background, residuals proj_model = self._represent_healpix_map( @@ -1199,15 +1160,11 @@ def display_fit(self, smoothing_kernel_sigma=0.1, display_colorbar=False): vmax = max(np.nanmax(proj_model), np.nanmax(proj_data)) # Plot model - images[0] = subs[i][0].imshow( - proj_model, origin="lower", vmin=vmin, vmax=vmax - ) + images[0] = subs[i][0].imshow(proj_model, origin="lower", vmin=vmin, vmax=vmax) subs[i][0].set_title("model, bin {}".format(data_analysis_bin.name)) # Plot data map - images[1] = subs[i][1].imshow( - proj_data, origin="lower", vmin=vmin, vmax=vmax - ) + images[1] = subs[i][1].imshow(proj_data, origin="lower", vmin=vmin, vmax=vmax) subs[i][1].set_title("excess, bin {}".format(data_analysis_bin.name)) # Plot background map. @@ -1330,9 +1287,7 @@ def _get_model_map(self, plane_id, n_pt_src, n_ext_src): raise ValueError(f"{plane_id} not a plane in the current model") model_map = SparseHealpix( - self._get_expectation( - self._maptree[plane_id], plane_id, n_pt_src, n_ext_src - ), + self._get_expectation(self._maptree[plane_id], plane_id, n_pt_src, n_ext_src), self._active_pixels[plane_id], self._maptree[plane_id].observation_map.nside, ) @@ -1406,17 +1361,13 @@ def _write_a_map(self, file_name, which, fluctuate=False, return_map=False): if return_map: return new_map_tree - def write_model_map( - self, file_name, poisson_fluctuate=False, test_return_map=False - ): + def write_model_map(self, file_name, poisson_fluctuate=False, test_return_map=False): """ This function writes the model map to a file. The interface is based off of HAWCLike for consistency """ if test_return_map: - log.warning( - "test_return_map=True should only be used for testing purposes!" - ) + log.warning("test_return_map=True should only be used for testing purposes!") return self._write_a_map(file_name, "model", poisson_fluctuate, test_return_map) def write_residual_map(self, file_name, test_return_map=False): @@ -1425,7 +1376,5 @@ def write_residual_map(self, file_name, test_return_map=False): The interface is based off of HAWCLike for consistency """ if test_return_map: - log.warning( - "test_return_map=True should only be used for testing purposes!" - ) + log.warning("test_return_map=True should only be used for testing purposes!") return self._write_a_map(file_name, "residual", False, test_return_map) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index 76f0e61..00bca4c 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/hawc_hal/maptree/from_hdf5_file.py @@ -1,19 +1,19 @@ from __future__ import absolute_import + import collections -from hawc_hal.serialize import Serialization from hawc_hal.region_of_interest import get_roi_from_dict - - +from hawc_hal.serialize import Serialization from threeML.io.logging import setup_logger + log = setup_logger(__name__) log.propagate = False -from ..healpix_handling import SparseHealpix, DenseHealpix +from ..healpix_handling import DenseHealpix, SparseHealpix from .data_analysis_bin import DataAnalysisBin -def from_hdf5_file(map_tree_file, roi): +def from_hdf5_file(map_tree_file, roi, transits): """ Create a MapTree object from a HDF5 file and a ROI. Do not use this directly, use map_tree_factory instead. @@ -25,12 +25,12 @@ def from_hdf5_file(map_tree_file, roi): # Read the data frames contained in the file with Serialization(map_tree_file) as serializer: - analysis_bins_df, _ = serializer.retrieve_pandas_object('/analysis_bins') - meta_df, _ = serializer.retrieve_pandas_object('/analysis_bins_meta') - roimap, roi_meta = serializer.retrieve_pandas_object('/ROI') - - if len(roimap)>0: - roi_meta['roimap']=roimap.values + analysis_bins_df, _ = serializer.retrieve_pandas_object("/analysis_bins") + meta_df, _ = serializer.retrieve_pandas_object("/analysis_bins_meta") + roimap, roi_meta = serializer.retrieve_pandas_object("/ROI") + + if len(roimap) > 0: + roi_meta["roimap"] = roimap.values # Let's see if the file contains the definition of an ROI if len(roi_meta) > 0: @@ -47,16 +47,19 @@ def from_hdf5_file(map_tree_file, roi): active_pixels_user = roi.active_pixels(1024) # This verifies that active_pixels_file is a superset (or equal) to the user-provided set - assert set(active_pixels_file) >= set(active_pixels_user), \ - "The ROI you provided (%s) is not a subset " \ + assert set(active_pixels_file) >= set(active_pixels_user), ( + "The ROI you provided (%s) is not a subset " "of the one contained in the file (%s)" % (roi, file_roi) + ) else: # The user has provided no ROI, but the file contains one. Let's issue a warning - log.warning("You did not provide any ROI but the map tree %s contains " - "only data within the ROI %s. " - "Only those will be used." % (map_tree_file, file_roi)) + log.warning( + "You did not provide any ROI but the map tree %s contains " + "only data within the ROI %s. " + "Only those will be used." % (map_tree_file, file_roi) + ) # Make a copy of the file ROI and use it as if the user provided that one roi = get_roi_from_dict(file_roi.to_dict()) @@ -69,6 +72,11 @@ def from_hdf5_file(map_tree_file, roi): data_analysis_bins = collections.OrderedDict() + if transits is not None: + n_transits = transits + else: + n_transits: float = meta_df["n_transits"].max() + for bin_name in bin_names: this_df = analysis_bins_df.loc[bin_name] @@ -77,31 +85,40 @@ def from_hdf5_file(map_tree_file, roi): if roi is not None: # Get the active pixels for this plane - active_pixels_user = roi.active_pixels(int(this_meta['nside'])) + active_pixels_user = roi.active_pixels(int(this_meta["nside"])) # Read only the pixels that the user wants - observation_hpx_map = SparseHealpix(this_df.loc[active_pixels_user, 'observation'].values, - active_pixels_user, this_meta['nside']) - background_hpx_map = SparseHealpix(this_df.loc[active_pixels_user, 'background'].values, - active_pixels_user, this_meta['nside']) + observation_hpx_map = SparseHealpix( + this_df.loc[active_pixels_user, "observation"].values, + active_pixels_user, + this_meta["nside"], + ) + background_hpx_map = SparseHealpix( + this_df.loc[active_pixels_user, "background"].values, + active_pixels_user, + this_meta["nside"], + ) else: # Full sky - observation_hpx_map = DenseHealpix(this_df.loc[:, 'observation'].values) - background_hpx_map = DenseHealpix(this_df.loc[:, 'background'].values) + observation_hpx_map = DenseHealpix(this_df.loc[:, "observation"].values) + background_hpx_map = DenseHealpix(this_df.loc[:, "background"].values) # 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'], - scheme='RING' if this_meta['scheme'] == 0 else 'NEST') + 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=n_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, n_transits diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index d70f10c..045127a 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/hawc_hal/maptree/from_root_file.py @@ -1,63 +1,32 @@ from __future__ import absolute_import -from builtins import map -from builtins import str -from builtins import range -from itertools import count + +import collections import os -from pathlib import Path import socket -import collections +from builtins import map, range, str +from itertools import count +from pathlib import Path from token import N_TOKENS + import healpy as hp -from matplotlib.style import library import numpy as np import uproot - +from matplotlib.style import library from threeML.io.file_utils import file_existing_and_readable, sanitize_filename - from threeML.io.logging import setup_logger log = setup_logger(__name__) log.propagate = False +from ..healpix_handling import DenseHealpix, SparseHealpix from ..region_of_interest import HealpixROIBase from .data_analysis_bin import DataAnalysisBin -from ..healpix_handling import SparseHealpix, DenseHealpix - -# ! As of right now, any code with ROOT functionality is commented, -# ! but should be removed in the future. -# def _get_bin_object(f, bin_name, suffix): - -# # kludge: some maptrees have bin labels in BinInfo like "0", "1", "2", but then -# # the nHit bin is actually "nHit00, nHit01, nHit02... others instead have -# # labels in BinInfo like "00", "01", "02", and still the nHit bin is nHit00, nHit01 -# # thus we need to add a 0 to the former case, but not add it to the latter case - -# # bin_label = "nHit0%s/%s" % (bin_name, suffix) -# bin_label = f"nHit0{bin_name}/{suffix}" - -# bin_tobject = f.Get(bin_label) - -# if not bin_tobject: - -# # Try the other way -# # bin_label = "nHit%s/%s" % (bin_name, suffix) -# bin_label = f"nHit{bin_name}/{suffix}" - -# bin_tobject = f.Get(bin_label) - -# if not bin_tobject: - -# # raise IOError("Could not read bin %s" % bin_label) -# raise IOError(f"Could not read bin {bin_label}") - -# return bin_tobject - def from_root_file( map_tree_file: Path, roi: HealpixROIBase, + transits: float, scheme: int = 0, ): """Create a Maptree object from a ROOT file and a ROI. @@ -112,14 +81,12 @@ def from_root_file( try: - # data_bins_labels = map_infile["BinInfo"]["name"].array().to_numpy() data_bins_labels = map_infile["BinInfo/name"].array().to_numpy() except ValueError: try: - # data_bins_labels = map_infile["BinInfo"]["id"].array().to_numpy() data_bins_labels = map_infile["BinInfo/id"].array().to_numpy() except ValueError as exc: @@ -135,12 +102,6 @@ def from_root_file( npix_cnt = map_infile[data_tree_prefix].array().to_numpy().size npix_bkg = map_infile[bkg_tree_prefix].array().to_numpy().size - # npix_cnt = ( - # map_infile[f"nHit{bin_name}"]["data"]["count"].array().to_numpy().size - # ) - # npix_bkg = ( - # map_infile[f"nHit{bin_name}"]["bkg"]["count"].array().to_numpy().size - # ) except uproot.KeyInFileError: data_tree_prefix = f"nHit0{bin_name}/data/count" @@ -149,21 +110,16 @@ def from_root_file( npix_cnt = map_infile[data_tree_prefix].array().to_numpy().size npix_bkg = map_infile[bkg_tree_prefix].array().to_numpy().size - # npix_cnt = ( - # map_infile[f"nHit0{bin_name}"]["data"]["count"].array().to_numpy().size - # ) - # npix_bkg = ( - # map_infile[f"nHit0{bin_name}"]["bkg"]["count"].array().to_numpy().size - # ) - # 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_durations = np.divide(map_infile["BinInfo"]["totalDuration"].array(), 24.0) maptree_durations = map_infile["BinInfo/totalDuration"].array() - n_durations = np.divide(maptree_durations, 24.0) - n_transits: float = max(n_durations) + if transits is not None: + n_transits = transits + else: + n_durations: np.ndarray = np.divide(maptree_durations, 24.0) + n_transits: float = max(n_durations) n_bins: int = data_bins_labels.shape[0] nside_cnt: int = hp.pixelfunc.npix2nside(npix_cnt) @@ -171,9 +127,7 @@ def from_root_file( # so far, a value of Nside of 1024 (perhaps will change) and a # RING HEALPix - assert ( - nside_cnt == nside_bkg - ), "Nside value needs to be the same for counts and bkg. maps" + assert nside_cnt == nside_bkg, "Nside value needs to be the same for counts and bkg. maps" assert scheme == 0, "NESTED HEALPix is not currently supported." @@ -184,9 +138,7 @@ def from_root_file( # HACK: simple way of reading the number of active pixels within # the define ROI if roi is not None: - active_pixels = roi.active_pixels( - nside_cnt, system="equatorial", ordering="RING" - ) + active_pixels = roi.active_pixels(nside_cnt, system="equatorial", ordering="RING") for pix_id in active_pixels: @@ -203,8 +155,6 @@ def from_root_file( counts = map_infile[data_tree_prefix].array().to_numpy() bkg = map_infile[bkg_tree_prefix].array().to_numpy() - # counts = map_infile[f"nHit{name}"]["data"]["count"].array().to_numpy() - # bkg = map_infile[f"nHit{name}"]["bkg"]["count"].array().to_numpy() except uproot.KeyInFileError: @@ -214,8 +164,6 @@ def from_root_file( counts = map_infile[data_tree_prefix].array().to_numpy() bkg = map_infile[bkg_tree_prefix].array().to_numpy() - # counts = map_infile[f"nHit0{name}"]["data"]["count"].array().to_numpy() - # bkg = map_infile[f"nHit0{name}"]["bkg"]["count"].array().to_numpy() # Read only elements within the ROI if roi is not None: @@ -226,9 +174,7 @@ def from_root_file( counts_hpx = SparseHealpix( counts[healpix_map_active > 0], active_pixels, nside_cnt ) - bkg_hpx = SparseHealpix( - bkg[healpix_map_active > 0], active_pixels, nside_cnt - ) + bkg_hpx = SparseHealpix(bkg[healpix_map_active > 0], active_pixels, nside_cnt) this_data_analysis_bin = DataAnalysisBin( name, @@ -242,21 +188,6 @@ def from_root_file( # Read the whole sky else: - # try: - - # counts = ( - # map_infile[f"nHit{name}"]["data"]["count"].array().to_numpy() - # ) - # bkg = map_infile[f"nHit{name}"]["bkg"]["count"].array().to_numpy() - - # except uproot.KeyInFileError: - - # # Sometimes, names of bins carry an extra zero - # counts = ( - # map_infile[f"nHit0{name}"]["data"]["count"].array().to_numpy() - # ) - # bkg = map_infile[f"nHit0{name}"]["bkg"]["count"].array().to_numpy() - this_data_analysis_bin = DataAnalysisBin( name, DenseHealpix(counts), @@ -268,165 +199,4 @@ def from_root_file( data_analysis_bins[name] = this_data_analysis_bin - # # Read map tree - # with open_ROOT_file(str(map_tree_file)) as f: - - # # Newer maps use "name" rather than "id" - # try: - - # data_bins_labels = list(root_numpy.tree2array(f.Get("BinInfo"), "name")) - - # except ValueError: - - # # Check to see if its an old style maptree - # try: - - # data_bins_labels = list(root_numpy.tree2array(f.Get("BinInfo"), "id")) - - # except ValueError: - - # # Give a useful error message - # raise ValueError("Maptree has no Branch: 'id' or 'name' ") - - # # If the old style, we need to make them strings - # 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) - - # n_bins = len(data_bins_labels) - - # # These are going to be Healpix maps, one for each data analysis bin_name - - # data_analysis_bins = collections.OrderedDict() - - # for i in range(n_bins): - - # name = data_bins_labels[i] - - # data_tobject = _get_bin_object(f, name, "data") - - # bkg_tobject = _get_bin_object(f, name, "bkg") - - # # Get ordering scheme - # nside = data_tobject.GetUserInfo().FindObject("Nside").GetVal() - # nside_bkg = bkg_tobject.GetUserInfo().FindObject("Nside").GetVal() - - # assert nside == nside_bkg - - # scheme = data_tobject.GetUserInfo().FindObject("Scheme").GetVal() - # scheme_bkg = bkg_tobject.GetUserInfo().FindObject("Scheme").GetVal() - - # assert scheme == scheme_bkg - - # assert scheme == 0, "NESTED scheme is not supported yet" - - # if roi is not None: - - # # Only read the elements in the ROI - - # active_pixels = roi.active_pixels( - # nside, system="equatorial", ordering="RING" - # ) - - # counts = _read_partial_tree(data_tobject, active_pixels) - # bkg = _read_partial_tree(bkg_tobject, active_pixels) - - # counts_hpx = SparseHealpix(counts, active_pixels, nside) - # bkg_hpx = SparseHealpix(bkg, active_pixels, nside) - - # this_data_analysis_bin = DataAnalysisBin( - # name, - # counts_hpx, - # bkg_hpx, - # active_pixels_ids=active_pixels, - # n_transits=n_transits, - # scheme="RING", - # ) - - # else: - - # # Read the entire sky. - - # counts = tree_to_ndarray(data_tobject, "count").astype(np.float64) - # bkg = tree_to_ndarray(bkg_tobject, "count").astype(np.float64) - - # this_data_analysis_bin = DataAnalysisBin( - # name, - # DenseHealpix(counts), - # DenseHealpix(bkg), - # active_pixels_ids=None, - # n_transits=n_transits, - # scheme="RING", - # ) - # - # data_analysis_bins[name] = this_data_analysis_bin - return data_analysis_bins - - -# def _read_partial_tree(ttree_instance, elements_to_read): - -# # Decide whether to use a smart loading scheme, or just loading the whole thing, based on the -# # number of elements to be read - -# from ..root_handler import ROOT, root_numpy, tree_to_ndarray - -# if elements_to_read.shape[0] < 500000: - -# # Use a smart loading scheme, where we read only the pixels we need - -# # The fastest method that I found is to create a TEventList, apply it to the -# # tree, get a copy of the subset and then use ttree2array - -# # Create TEventList -# entrylist = ROOT.TEntryList() - -# # Add the selections -# _ = list(map(entrylist.Enter, elements_to_read)) - -# # Apply the EntryList to the tree -# ttree_instance.SetEntryList(entrylist) - -# # Get copy of the subset -# # We need to create a dumb TFile to silence a lot of warnings from ROOT -# # Get a filename for this process -# dumb_tfile_name = "__dumb_tfile_%s_%s.root" % ( -# os.getpid(), -# socket.gethostname(), -# ) -# dumb_tfile = ROOT.TFile(dumb_tfile_name, "RECREATE") -# new_tree = ttree_instance.CopyTree("") - -# # Actually read it from disk -# partial_map = root_numpy.tree2array(new_tree, "count").astype(np.float64) - -# # Now reset the entry list -# ttree_instance.SetEntryList(0) - -# dumb_tfile.Close() -# os.remove(dumb_tfile_name) - -# else: - -# # The smart scheme starts to become slower than the brute force approach, so let's read the whole thing -# partial_map = tree_to_ndarray(ttree_instance, "count").astype(np.float64) - -# assert ( -# partial_map.shape[0] >= elements_to_read.shape[0] -# ), "Trying to read more pixels than present in TTree" - -# # Unless we have read the whole sky, let's remove the pixels we shouldn't have read - -# if elements_to_read.shape[0] != partial_map.shape[0]: - -# # Now let's remove the pixels we shouldn't have read -# partial_map = partial_map[elements_to_read] - -# return partial_map.astype(np.float64) + return data_analysis_bins, n_transits diff --git a/hawc_hal/maptree/map_tree.py b/hawc_hal/maptree/map_tree.py index 35b952a..071857f 100644 --- a/hawc_hal/maptree/map_tree.py +++ b/hawc_hal/maptree/map_tree.py @@ -1,55 +1,67 @@ -from __future__ import division -from __future__ import absolute_import -from builtins import object -from past.utils import old_div +from __future__ import absolute_import, division + import os -import uproot +from builtins import object + import numpy as np import pandas as pd - -from threeML.io.rich_display import display +import uproot +from past.utils import old_div from threeML.io.file_utils import sanitize_filename from threeML.io.logging import setup_logger +from threeML.io.rich_display import display log = setup_logger(__name__) log.propagate = False +import astropy.units as u + from ..serialize import Serialization -from .from_root_file import from_root_file from .from_hdf5_file import from_hdf5_file - -import astropy.units as u +from .from_root_file import from_root_file -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 is not None: + self._n_transits = n_transits + else: + # find the transits in the map + transits = [] + for analysis_bin_key in self._analysis_bins: + transits.append(self._analysis_bins[analysis_bin_key].n_transits) + + # get the largest value + 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): + def from_hdf5(cls, map_tree_file, roi, n_transits): - data_analysis_bins = from_hdf5_file(map_tree_file, roi) + data_analysis_bins, transits = from_hdf5_file(map_tree_file, roi, n_transits) - return cls(data_analysis_bins, roi) + 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. @@ -58,9 +70,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): """ @@ -99,6 +111,10 @@ def __len__(self): return len(self._analysis_bins) + @property + def n_transits(self): + return self._n_transits + @property def analysis_bins_labels(self): @@ -109,12 +125,8 @@ def display(self): df = pd.DataFrame() df["Bin"] = list(self._analysis_bins.keys()) - df["Nside"] = [ - self._analysis_bins[bin_id].nside for bin_id in self._analysis_bins - ] - df["Scheme"] = [ - self._analysis_bins[bin_id].scheme for bin_id in self._analysis_bins - ] + df["Nside"] = [self._analysis_bins[bin_id].nside for bin_id in self._analysis_bins] + df["Scheme"] = [self._analysis_bins[bin_id].scheme for bin_id in self._analysis_bins] # Compute observed counts, background counts, how many pixels we have in the ROI and # the sky area they cover @@ -181,9 +193,9 @@ def write(self, filename): analysis_bin = self._analysis_bins[bin_id] - assert ( - bin_id == analysis_bin.name - ), "Bin name inconsistency: {} != {}".format(bin_id, analysis_bin.name) + assert bin_id == analysis_bin.name, "Bin name inconsistency: {} != {}".format( + bin_id, analysis_bin.name + ) multi_index_keys.append(analysis_bin.name) @@ -211,9 +223,7 @@ def write(self, filename): else: - serializer.store_pandas_object( - "/ROI", pd.Series(), **self._roi.to_dict() - ) + serializer.store_pandas_object("/ROI", pd.Series(), **self._roi.to_dict()) else: diff --git a/hawc_hal/psf_fast/psf_wrapper.py b/hawc_hal/psf_fast/psf_wrapper.py index c65035d..1279e5d 100644 --- a/hawc_hal/psf_fast/psf_wrapper.py +++ b/hawc_hal/psf_fast/psf_wrapper.py @@ -1,13 +1,14 @@ from __future__ import division -from builtins import zip -from builtins import object -from past.utils import old_div + import copy +from builtins import object, zip + import numpy as np +import pandas as pd import scipy.integrate import scipy.interpolate import scipy.optimize -import pandas as pd +from past.utils import old_div _INTEGRAL_OUTER_RADIUS = 15.0 @@ -29,9 +30,7 @@ def __init__(self, xs, ys, brightness_interp_x=None, brightness_interp_y=None): # Memorize the total integral (will use it for normalization) - self._total_integral = self._psf_interpolated.integral( - self._xs[0], _INTEGRAL_OUTER_RADIUS - ) + self._total_integral = self._psf_interpolated.integral(self._xs[0], _INTEGRAL_OUTER_RADIUS) # Now compute the truncation radius, which is a very conservative measurement # of the size of the PSF @@ -91,9 +90,7 @@ def find_eef_radius(self, fraction): f = lambda r: fraction - old_div(self.integral(1e-4, r), self._total_integral) - radius, status = scipy.optimize.brentq( - f, 0.005, _INTEGRAL_OUTER_RADIUS, full_output=True - ) + radius, status = scipy.optimize.brentq(f, 0.005, _INTEGRAL_OUTER_RADIUS, full_output=True) assert status.converged, "Brentq did not converged" @@ -134,9 +131,7 @@ def combine_with_other_psf(self, other_psf, w1, w2): new_ys = w1 * self.ys + w2 * other_psf.ys # Also weight the brightness interpolation points - new_br_interp_y = ( - w1 * self._brightness_interp_y + w2 * other_psf._brightness_interp_y - ) + new_br_interp_y = w1 * self._brightness_interp_y + w2 * other_psf._brightness_interp_y return PSFWrapper( self.xs, @@ -228,29 +223,6 @@ def psf_function(ang_dist: float): return new_instance - # @classmethod - # def from_TF1(cls, tf1_instance): - # - # # Annoyingly, some PSFs for some Dec bins (at large Zenith angles) have - # # zero or negative integrals, i.e., they are useless. Return an unusable PSF - # # object in that case - # if tf1_instance.Integral(0, _INTEGRAL_OUTER_RADIUS) <= 0.0: - # - # return InvalidPSF() - # - # # Make interpolation - # xs = np.logspace(-3, np.log10(_INTEGRAL_OUTER_RADIUS), 500) - # ys = np.array([tf1_instance.Eval(x) for x in xs], float) - # - # assert np.all(np.isfinite(xs)) - # assert np.all(np.isfinite(ys)) - # - # instance = cls(xs, ys) - # - # instance._tf1 = tf1_instance.Clone() - # - # return instance - def integral(self, a, b): return self._psf_interpolated.integral(a, b) diff --git a/hawc_hal/response/response.py b/hawc_hal/response/response.py index 1999f69..7f37bff 100644 --- a/hawc_hal/response/response.py +++ b/hawc_hal/response/response.py @@ -1,23 +1,21 @@ -from __future__ import division -from __future__ import absolute_import -from builtins import zip -from builtins import range -from builtins import object +from __future__ import absolute_import, division + +import collections +import os +from builtins import object, range, zip from multiprocessing.managers import ValueProxy from pathlib import Path -from past.utils import old_div -import uproot + import healpy as hp import numpy as np import pandas as pd -import os -import collections +import uproot +from past.utils import old_div +from threeML.io.file_utils import file_existing_and_readable, sanitize_filename +from threeML.io.logging import setup_logger from ..serialize import Serialization -from threeML.io.file_utils import file_existing_and_readable, sanitize_filename - -from threeML.io.logging import setup_logger log = setup_logger(__name__) log.propagate = False @@ -53,7 +51,7 @@ def hawc_response_factory(response_file_name): new_instance = HAWCResponse.from_root_file(response_file_name) - elif extension in ['.hd5', '.hdf5', '.hdf']: + elif extension in [".hd5", ".hdf5", ".hdf"]: new_instance = HAWCResponse.from_hdf5(response_file_name) @@ -63,7 +61,7 @@ def hawc_response_factory(response_file_name): f"Extension {extension} for response file {response_file_name} not recognized." ) # raise NotImplementedError("Extension %s for response file %s not recognized." % (extension, - # response_file_name)) + # response_file_name)) _instances[response_file_name] = new_instance @@ -73,13 +71,12 @@ def hawc_response_factory(response_file_name): class HAWCResponse(object): - def __init__(self, response_file_name, dec_bins, response_bins): self._response_file_name = response_file_name self._dec_bins = dec_bins self._response_bins = response_bins - + if len(dec_bins) < 2: # log.warning("Only {0} dec bins given in {1}, will not try to interpolate.".format(len(dec_bins), response_file_name)) log.warning( @@ -101,11 +98,11 @@ def from_hdf5(cls, response_file_name): response_bins = collections.OrderedDict() - with Serialization(response_file_name, mode='r') as serializer: + with Serialization(response_file_name, mode="r") as serializer: - meta_dfs, _ = serializer.retrieve_pandas_object('/dec_bins_definition') - effarea_dfs, _ = serializer.retrieve_pandas_object('/effective_area') - psf_dfs, _ = serializer.retrieve_pandas_object('/psf') + meta_dfs, _ = serializer.retrieve_pandas_object("/dec_bins_definition") + effarea_dfs, _ = serializer.retrieve_pandas_object("/effective_area") + psf_dfs, _ = serializer.retrieve_pandas_object("/psf") declination_centers = effarea_dfs.index.levels[0] energy_bins = effarea_dfs.index.levels[1] @@ -121,9 +118,9 @@ def from_hdf5(cls, response_file_name): these_meta = meta_dfs.loc[dec_center, energy_bin] - min_dec = these_meta['min_dec'] - max_dec = these_meta['max_dec'] - dec_center_ = these_meta['declination_center'] + min_dec = these_meta["min_dec"] + max_dec = these_meta["max_dec"] + dec_center_ = these_meta["declination_center"] assert dec_center_ == dec_center, "Response is corrupted" @@ -141,27 +138,37 @@ def from_hdf5(cls, response_file_name): assert min_dec == min_decs[-1], "Response is corrupted" assert max_dec == max_decs[-1], "Response is corrupted" - sim_n_sig_events = these_meta['n_sim_signal_events'] - sim_n_bg_events = these_meta['n_sim_bkg_events'] + sim_n_sig_events = these_meta["n_sim_signal_events"] + sim_n_bg_events = these_meta["n_sim_bkg_events"] this_effarea_df = effarea_dfs.loc[dec_center, energy_bin] - sim_energy_bin_low = this_effarea_df.loc[:, 'sim_energy_bin_low'].values - sim_energy_bin_centers = this_effarea_df.loc[:, 'sim_energy_bin_centers'].values - sim_energy_bin_hi = this_effarea_df.loc[:, 'sim_energy_bin_hi'].values - sim_differential_photon_fluxes = this_effarea_df.loc[:, 'sim_differential_photon_fluxes'].values - sim_signal_events_per_bin = this_effarea_df.loc[:, 'sim_signal_events_per_bin'].values + sim_energy_bin_low = this_effarea_df.loc[:, "sim_energy_bin_low"].values + sim_energy_bin_centers = this_effarea_df.loc[:, "sim_energy_bin_centers"].values + sim_energy_bin_hi = this_effarea_df.loc[:, "sim_energy_bin_hi"].values + sim_differential_photon_fluxes = this_effarea_df.loc[ + :, "sim_differential_photon_fluxes" + ].values + sim_signal_events_per_bin = this_effarea_df.loc[ + :, "sim_signal_events_per_bin" + ].values this_psf = PSFWrapper.from_pandas(psf_dfs.loc[dec_center, energy_bin, :]) - this_response_bin = ResponseBin(energy_bin, min_dec, max_dec, dec_center, - sim_n_sig_events, sim_n_bg_events, - sim_energy_bin_low, - sim_energy_bin_centers, - sim_energy_bin_hi, - sim_differential_photon_fluxes, - sim_signal_events_per_bin, - this_psf) + this_response_bin = ResponseBin( + energy_bin, + min_dec, + max_dec, + dec_center, + sim_n_sig_events, + sim_n_bg_events, + sim_energy_bin_low, + sim_energy_bin_centers, + sim_energy_bin_hi, + sim_differential_photon_fluxes, + sim_signal_events_per_bin, + this_psf, + ) these_response_bins[energy_bin] = this_response_bin @@ -208,10 +215,7 @@ def from_root_file(cls, response_file_name: Path): # NOTE: The LogLogSpectrum parameters are extracted from the response file # There is no way to evaluate the loglogspectrum function with uproot, so the # parameters are passed for evaluation in response_bin.py - # 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() - + log_log_params = response_file["LogLogSpectrum"].member("fParams") log_log_shape = response_file["LogLogSpectrum"].member("fTitle") dec_bins_lower_edge = response_file["DecBins/lowerEdge"].array().to_numpy() @@ -222,14 +226,12 @@ def from_root_file(cls, response_file_name: Path): try: - # response_bin_ids = response_file["AnalysisBins"]["name"].array().to_numpy() response_bin_ids = response_file["AnalysisBins/name"].array().to_numpy() except uproot.KeyInFileError: try: - # response_bin_ids = response_file["AnalysisBins"]["id"].array().to_numpy() response_bin_ids = response_file["AnalysisBins/id"].array().to_numpy() except uproot.KeyInFileError: @@ -245,7 +247,7 @@ 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)): this_response_bins = collections.OrderedDict() @@ -261,7 +263,7 @@ def from_root_file(cls, response_file_name: Path): response_bins_ids = list(range(n_energy_bins)) for response_bin_id in response_bin_ids: - + this_response_bin = ResponseBin.from_ttree( response_file, dec_id, @@ -270,7 +272,8 @@ def from_root_file(cls, response_file_name: Path): log_log_shape, min_dec, dec_center, - max_dec) + max_dec, + ) this_response_bins[response_bin_id] = this_response_bin @@ -278,92 +281,15 @@ def from_root_file(cls, response_file_name: Path): del response_file - # with open_ROOT_file( str(response_file_name) ) as root_file: - - # # Get the name of the trees - # object_names = get_list_of_keys(root_file) - - # # Make sure we have all the things we need - - # assert 'LogLogSpectrum' in object_names - # assert 'DecBins' in object_names - # assert 'AnalysisBins' in object_names - - # # Read spectrum used during the simulation - # log_log_spectrum = root_file.Get("LogLogSpectrum") - - # # Get the analysis bins definition - # dec_bins_ = tree_to_ndarray(root_file.Get("DecBins")) - - # dec_bins_lower_edge = dec_bins_['lowerEdge'] # type: np.ndarray - # dec_bins_upper_edge = dec_bins_['upperEdge'] # type: np.ndarray - # dec_bins_center = dec_bins_['simdec'] # type: np.ndarray - - # dec_bins = list(zip(dec_bins_lower_edge, dec_bins_center, dec_bins_upper_edge)) - - # # Read in the ids of the response bins ("analysis bins" in LiFF jargon) - # try: - - # response_bins_ids = tree_to_ndarray(root_file.Get("AnalysisBins"), "name") # type: np.ndarray - - # except ValueError: - - # try: - - # response_bins_ids = tree_to_ndarray(root_file.Get("AnalysisBins"), "id") # type: np.ndarray - - # except ValueError: - - # # Some old response files (or energy responses) have no "name" branch - # log.warning("Response %s has no AnalysisBins 'id' or 'name' branch. " - # "Will try with default names" % response_file_name) - - # response_bins_ids = None - - # response_bins_ids = response_bins_ids.astype(str) - - # # 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)): - - # this_response_bins = collections.OrderedDict() - - # min_dec, dec_center, max_dec = dec_bins[dec_id] - - # # If we couldn't get the reponse_bins_ids above, let's use the default names - # if response_bins_ids is None: - - # # Default are just integers. let's read how many nHit bins are from the first dec bin - # dec_id_label = "dec_%02i" % dec_id - - # n_energy_bins = root_file.Get(dec_id_label).GetNkeys() - - # response_bins_ids = list(range(n_energy_bins)) - - # for response_bin_id in response_bins_ids: - # this_response_bin = ResponseBin.from_ttree(root_file, dec_id, response_bin_id, log_log_spectrum, - # min_dec, dec_center, max_dec) - - # this_response_bins[response_bin_id] = this_response_bin - - # response_bins[dec_bins[dec_id][1]] = this_response_bins - - # # Now the file is closed. Let's explicitly remove f so we are sure it is freed - # del root_file - - # Instance the class and return it - # return instance return cls(response_file_name, dec_bins, response_bins) - def get_response_dec_bin(self, dec, interpolate=False): """Get the response for the provided declination bin, optionally interpolating the PSF Args: dec (float): Declination where the response is desired - interpolate (bool, optional): If True, PSF is interpolated between the two + interpolate (bool, optional): If True, PSF is interpolated between the two closest response bins. Defaults to False. Returns: @@ -374,7 +300,7 @@ def get_response_dec_bin(self, dec, interpolate=False): dec_bins_keys = list(self._response_bins.keys()) dec_bins_by_distance = sorted(dec_bins_keys, key=lambda x: abs(x - dec)) - #never try to interpolate if only one dec bin is given + # never try to interpolate if only one dec bin is given if len(dec_bins_keys) < 2: interpolate = False @@ -410,7 +336,9 @@ def get_response_dec_bin(self, dec, interpolate=False): new_responses = collections.OrderedDict() for bin_id in energy_bins_one: - this_new_response = energy_bins_one[bin_id].combine_with_weights(energy_bins_two[bin_id], dec, w1, w2) + this_new_response = energy_bins_one[bin_id].combine_with_weights( + energy_bins_two[bin_id], dec, w1, w2 + ) new_responses[bin_id] = this_new_response @@ -481,9 +409,10 @@ def write(self, filename): effarea_dfs.append(this_effarea_df) psf_dfs.append(this_psf_df) # assert bin_id == response_bin.name, \ - # 'Bin name inconsistency: {} != {}'.format(bin_id, response_bin.name) - assert bin_id == response_bin.name, \ - f'Bin name inconsistency: {bin_id} != {response_bin.name}' + # 'Bin name inconsistency: {} != {}'.format(bin_id, response_bin.name) + assert ( + bin_id == response_bin.name + ), f"Bin name inconsistency: {bin_id} != {response_bin.name}" multi_index_keys.append((dec_center, response_bin.name)) all_metas.append(pd.Series(this_meta)) @@ -493,8 +422,8 @@ def write(self, filename): meta_df = pd.concat(all_metas, axis=1, keys=multi_index_keys).T # Now write the 4 dataframes to file - with Serialization(filename, mode='w') as serializer: + with Serialization(filename, mode="w") as serializer: - serializer.store_pandas_object('/dec_bins_definition', meta_df) - serializer.store_pandas_object('/effective_area', effarea_df) - serializer.store_pandas_object('/psf', psf_df) + serializer.store_pandas_object("/dec_bins_definition", meta_df) + serializer.store_pandas_object("/effective_area", effarea_df) + serializer.store_pandas_object("/psf", psf_df) diff --git a/hawc_hal/response/response_bin.py b/hawc_hal/response/response_bin.py index 9d64598..06297f6 100644 --- a/hawc_hal/response/response_bin.py +++ b/hawc_hal/response/response_bin.py @@ -50,29 +50,20 @@ def _get_en_th1d( prefix: str, ): - # en_sig_label = "En%s_dec%i_nh%s" % (prefix, dec_id, analysis_bin_id) en_sig_label = f"En{prefix}_dec{dec_id}_nh{analysis_bin_id}" - # this_en_th1d = open_ttree.FindObjectAny(en_sig_label) try: hist_prefix = f"dec_{dec_id:02d}/nh_{analysis_bin_id}/{en_sig_label}" this_en_th1d = open_ttree[hist_prefix].to_hist() - # this_en_th1d = open_ttree[f"dec_{dec_id:02d}"][f"nh_{analysis_bin_id}"][ - # en_sig_label - # ].to_hist() except uproot.KeyInFileError: hist_prefix = f"dec_{dec_id:02d}/nh_0{analysis_bin_id}/{en_sig_label}" this_en_th1d = open_ttree[hist_prefix].to_hist() - # this_en_th1d = open_ttree[f"dec_{dec_id:02d}"][f"nh_0{analysis_bin_id}"][ - # en_sig_label - # ].to_hist() if not this_en_th1d: raise IOError(f"Could not find TH1D named {en_sig_label}.") - # raise IOError("Could not find TH1D named %s." % en_sig_label) return this_en_th1d @@ -121,20 +112,15 @@ def log_log_spectrum(log_energy: float, parameters: np.ndarray): return ( np.log10(parameters[0]) - parameters[1] * log_energy - - np.log10(np.exp(1.0)) - * np.power(10.0, log_energy - np.log10(parameters[2])) + - np.log10(np.exp(1.0)) * np.power(10.0, log_energy - np.log10(parameters[2])) ) else: raise ValueError("Unknown spectral shape.") - this_en_sig_th1d = cls._get_en_th1d( - root_file, dec_id, analysis_bin_id, prefix="Sig" - ) + this_en_sig_th1d = cls._get_en_th1d(root_file, dec_id, analysis_bin_id, prefix="Sig") - this_en_bg_th1d = cls._get_en_th1d( - root_file, dec_id, analysis_bin_id, prefix="Bg" - ) + this_en_bg_th1d = cls._get_en_th1d(root_file, dec_id, analysis_bin_id, prefix="Bg") total_bins = this_en_sig_th1d.shape[0] sim_energy_bin_low = np.zeros(total_bins) @@ -177,22 +163,11 @@ def log_log_spectrum(log_energy: float, parameters: np.ndarray): try: psf_prefix = f"dec_{dec_id:02d}/nh_{analysis_bin_id}" - psf_tf1_metadata = root_file[ - f"{psf_prefix}/PSF_dec{dec_id}_nh{analysis_bin_id}_fit" - ] - # psf_tf1_metadata = root_file[f"dec_{dec_id:02d}"][f"nh_{analysis_bin_id}"][ - # f"PSF_dec{dec_id}_nh{analysis_bin_id}_fit" - # ] + psf_tf1_metadata = root_file[f"{psf_prefix}/PSF_dec{dec_id}_nh{analysis_bin_id}_fit"] except uproot.KeyInFileError: psf_prefix = f"dec_{dec_id:02d}/nh_0{analysis_bin_id}" - psf_tf1_metadata = root_file[ - f"{psf_prefix}/PSF_dec{dec_id}_nh{analysis_bin_id}_fit" - ] - - # psf_tf1_metadata = root_file[f"dec_{dec_id:02d}"][f"nh_0{analysis_bin_id}"][ - # f"PSF_dec{dec_id}_nh{analysis_bin_id}_fit" - # ] + psf_tf1_metadata = root_file[f"{psf_prefix}/PSF_dec{dec_id}_nh{analysis_bin_id}_fit"] psf_tf1_fparams = psf_tf1_metadata.member("fParams") @@ -213,94 +188,6 @@ def log_log_spectrum(log_energy: float, parameters: np.ndarray): psf_fun, ) - # @classmethod - # def from_ttree( - # cls, - # open_ttree, - # dec_id, - # analysis_bin_id, - # log_log_spectrum, - # min_dec, - # dec_center, - # max_dec, - # ): - - # from ..root_handler import ROOT - - # # Read the histogram of the simulated events detected in this bin_name - # # NOTE: we do not copy this TH1D instance because we won't use it after the - # # file is closed - - # this_en_sig_th1d = cls._get_en_th1d(open_ttree, dec_id, analysis_bin_id, "Sig") - - # # The sum of the histogram is the total number of simulated events detected - # # in this analysis bin_name - # sim_n_sig_events = this_en_sig_th1d.Integral() - - # # Now let's see what has been simulated, i.e., the differential flux - # # at the center of each bin_name of the en_sig histogram - # sim_energy_bin_low = np.zeros(this_en_sig_th1d.GetNbinsX()) - # sim_energy_bin_centers = np.zeros(this_en_sig_th1d.GetNbinsX()) - # sim_energy_bin_hi = np.zeros(this_en_sig_th1d.GetNbinsX()) - # sim_signal_events_per_bin = np.zeros_like(sim_energy_bin_centers) - # sim_differential_photon_fluxes = np.zeros_like(sim_energy_bin_centers) - - # for i in range(sim_energy_bin_centers.shape[0]): - # # Remember: bin_name 0 is the underflow bin_name, that is why there - # # is a "i+1" and not just "i" - # bin_lo = this_en_sig_th1d.GetBinLowEdge(i + 1) - # bin_center = this_en_sig_th1d.GetBinCenter(i + 1) - # bin_hi = this_en_sig_th1d.GetBinWidth(i + 1) + bin_lo - - # # Store the center of the logarithmic bin_name - # sim_energy_bin_low[i] = 10**bin_lo # TeV - # sim_energy_bin_centers[i] = 10**bin_center # TeV - # sim_energy_bin_hi[i] = 10**bin_hi # TeV - - # # Get from the simulated spectrum the value of the differential flux - # # at the center energy - # sim_differential_photon_fluxes[i] = 10 ** log_log_spectrum.Eval( - # bin_center - # ) # TeV^-1 cm^-1 s^-1 - - # # Get from the histogram the detected events in each log-energy bin_name - # sim_signal_events_per_bin[i] = this_en_sig_th1d.GetBinContent(i + 1) - - # # Read the histogram of the bkg events detected in this bin_name - # # NOTE: we do not copy this TH1D instance because we won't use it after the - # # file is closed - - # this_en_bg_th1d = cls._get_en_th1d(open_ttree, dec_id, analysis_bin_id, "Bg") - - # # The sum of the histogram is the total number of simulated events detected - # # in this analysis bin_name - # sim_n_bg_events = this_en_bg_th1d.Integral() - - # # Now read the various TF1(s) for PSF, signal and background - - # # Read the PSF and make a copy (so it will stay when we close the file) - - # psf_label_tf1 = "PSF_dec%i_nh%s_fit" % (dec_id, analysis_bin_id) - - # tf1 = open_ttree.FindObjectAny(psf_label_tf1) - - # psf_fun = PSFWrapper.from_TF1(tf1) - - # return cls( - # analysis_bin_id, - # min_dec, - # max_dec, - # dec_center, - # sim_n_sig_events, - # sim_n_bg_events, - # sim_energy_bin_low, - # sim_energy_bin_centers, - # sim_energy_bin_hi, - # sim_differential_photon_fluxes, - # sim_signal_events_per_bin, - # psf_fun, - # ) - def to_pandas(self): """Save the information from Response file into a pandas.DataFrame @@ -354,9 +241,7 @@ def combine_with_weights(self, other_response_bin, dec_center, w1, w2): n_sim_signal_events = ( w1 * self._sim_n_sig_events + w2 * other_response_bin._sim_n_sig_events ) - n_sim_bkg_events = ( - w1 * self._sim_n_bg_events + w2 * other_response_bin._sim_n_bg_events - ) + n_sim_bkg_events = w1 * self._sim_n_bg_events + w2 * other_response_bin._sim_n_bg_events # We assume that the bin centers are the same assert np.allclose( diff --git a/hawc_hal/tests/conftest.py b/hawc_hal/tests/conftest.py index ec21a3e..f3706a0 100644 --- a/hawc_hal/tests/conftest.py +++ b/hawc_hal/tests/conftest.py @@ -1,49 +1,40 @@ -from __future__ import division -from __future__ import print_function -from past.utils import old_div +from __future__ import division, print_function + import matplotlib +from past.utils import old_div + matplotlib.use("Agg") import os -os.environ['OMP_NUM_THREADS'] = '1' -os.environ['MKL_NUM_THREADS'] = '1' -os.environ['NUMEXPR_NUM_THREADS'] = '1' +os.environ["OMP_NUM_THREADS"] = "1" +os.environ["MKL_NUM_THREADS"] = "1" +os.environ["NUMEXPR_NUM_THREADS"] = "1" -import pytest -from hawc_hal import HealpixConeROI -from threeML import * -from hawc_hal import HAL import time + import numpy as np +import pytest +from hawc_hal import HAL, HealpixConeROI +from threeML import * # Get data path -test_data_path = os.path.join(os.path.abspath(os.path.dirname(__file__)), 'data') +test_data_path = os.path.join(os.path.abspath(os.path.dirname(__file__)), "data") -def fit_point_source(roi, - maptree, - response, - point_source_model, - liff=False): +def fit_point_source(roi, maptree, response, point_source_model, liff=False): data_radius = roi.data_radius.to("deg").value if not liff: # This is a 3ML plugin - hawc = HAL("HAWC", - maptree, - response, - roi) + hawc = HAL("HAWC", maptree, response, roi) else: from threeML import HAWCLike - hawc = HAWCLike("HAWC", - maptree, - response, - fullsky=True) + hawc = HAWCLike("HAWC", maptree, response, fullsky=True) ra_roi, dec_roi = roi.ra_dec_center @@ -51,7 +42,8 @@ def fit_point_source(roi, hawc.set_active_measurements(1, 9) - if not liff: hawc.display() + if not liff: + hawc.display() data = DataList(hawc) @@ -73,7 +65,6 @@ def fit_point_source(roi, hawc.display_fit(display_colorbar=True).savefig("display_fit.png") - return param_df, like_df @@ -93,6 +84,18 @@ def check_map_trees(m1, m2): assert p1.n_transits == p2.n_transits +def check_n_transits(maptree, transits): + # Should give maptree object and expected transits + + # check maptree is assigned proper transits + assert maptree.n_transits == transits + + # check each bin as well + for bin_key in maptree: + the_bin = maptree[bin_key] + assert the_bin.n_transits == transits + + def check_responses(r1, r2): assert len(r1.response_bins) == len(r2.response_bins) @@ -121,7 +124,9 @@ def check_responses(r1, r2): assert np.allclose(rb1.sim_energy_bin_centers, rb2.sim_energy_bin_centers) - assert np.allclose(rb1.sim_differential_photon_fluxes, rb2.sim_differential_photon_fluxes) + assert np.allclose( + rb1.sim_differential_photon_fluxes, rb2.sim_differential_photon_fluxes + ) assert np.allclose(rb1.sim_signal_events_per_bin, rb2.sim_signal_events_per_bin) @@ -149,10 +154,12 @@ def roi(): # NOTE: the center of the ROI is not exactly on the source. This is on purpose, to make sure that we are # doing the model map with the right orientation - roi = HealpixConeROI(data_radius=data_radius, - model_radius=model_radius, - ra=ra_sim_source - 1.0, - dec=dec_sim_source + 0.5) + roi = HealpixConeROI( + data_radius=data_radius, + model_radius=model_radius, + ra=ra_sim_source - 1.0, + dec=dec_sim_source + 0.5, + ) return roi @@ -160,15 +167,12 @@ def roi(): @pytest.fixture(scope="session", autouse=True) def geminga_roi(): - ra_c, dec_c, rad = 101.7, 16, 9. + ra_c, dec_c, rad = 101.7, 16, 9.0 # NOTE: the center of the ROI is not exactly on the source. This is on purpose, to make sure that we are # doing the model map with the right orientation - roi = HealpixConeROI(data_radius=rad, - model_radius=rad + 15.0, - ra=ra_c, - dec=dec_c) + roi = HealpixConeROI(data_radius=rad, model_radius=rad + 15.0, ra=ra_c, dec=dec_c) return roi @@ -176,25 +180,25 @@ def geminga_roi(): @pytest.fixture(scope="session", autouse=True) def geminga_maptree(): - return os.path.join(test_data_path, 'geminga_maptree.root') + return os.path.join(test_data_path, "geminga_maptree.root") @pytest.fixture(scope="session", autouse=True) def geminga_response(): - return os.path.join(test_data_path, 'geminga_response.root') + return os.path.join(test_data_path, "geminga_response.root") @pytest.fixture(scope="session", autouse=True) def maptree(): - return os.path.join(test_data_path, 'zebra_simulated_source_mt_roi.hd5') + return os.path.join(test_data_path, "zebra_simulated_source_mt_roi.hd5") @pytest.fixture(scope="session", autouse=True) def response(): - return os.path.join(test_data_path, 'detector_response.hd5') + return os.path.join(test_data_path, "detector_response.hd5") @pytest.fixture(scope="module", autouse=True) @@ -217,7 +221,7 @@ def point_source_model(ra=83.6279, dec=22.14): spectrum.piv.fix = True - spectrum.K = old_div(3.15e-11, (u.TeV * u.cm ** 2 * u.s)) # norm (in 1/(keV cm2 s)) + spectrum.K = old_div(3.15e-11, (u.TeV * u.cm**2 * u.s)) # norm (in 1/(keV cm2 s)) spectrum.K.bounds = (1e-23, 1e-17) # without units energies are in keV spectrum.index = -2.0 @@ -230,4 +234,3 @@ def point_source_model(ra=83.6279, dec=22.14): model = Model(source) return model - diff --git a/hawc_hal/tests/test_n_transits.py b/hawc_hal/tests/test_n_transits.py new file mode 100644 index 0000000..7a26400 --- /dev/null +++ b/hawc_hal/tests/test_n_transits.py @@ -0,0 +1,33 @@ +import os + +import pytest +from hawc_hal import HealpixConeROI +from hawc_hal.maptree.map_tree import map_tree_factory + +from conftest import check_n_transits + + +def test_transits(maptree, roi): + """Specify the number of transits to use with maptree + + Args: + maptree (path.Path): Maptree in either ROOT or HDF5 format + roi (roi): Region of interest for analysis + """ + roi_ = roi + + # Case 1: specify number of transits + + n_transits = 777.7 + maptree_ntransits = map_tree_factory(maptree,roi_,n_transits) + + # does the maptree return the specified transits? + check_n_transits(maptree_ntransits, n_transits) + + # Case 2: number of transits is not specified (use value from maptree) + maptree_unspecifed = map_tree_factory(maptree, roi_) + n_transits = maptree_unspecifed.n_transits + + check_n_transits(maptree_unspecifed, n_transits) + + \ No newline at end of file From 568e132f42488182b5a479dffb7c7a4b6d912b79 Mon Sep 17 00:00:00 2001 From: Ramiro Torres-Escobedo Date: Thu, 15 Sep 2022 17:39:11 -0500 Subject: [PATCH 2/6] scale counts in data & bkg for transit number --- hawc_hal/maptree/from_root_file.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index 045127a..b6b90b5 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/hawc_hal/maptree/from_root_file.py @@ -115,10 +115,14 @@ def from_root_file( # estimate of the livetime is the maximum of n_transits, which normally # happen in the bins with high statistic maptree_durations = map_infile["BinInfo/totalDuration"].array() + n_durations: np.ndarray = np.divide(maptree_durations, 24.0) + if transits is not None: + n_transits = transits + else: - n_durations: np.ndarray = np.divide(maptree_durations, 24.0) + n_transits: float = max(n_durations) n_bins: int = data_bins_labels.shape[0] @@ -138,6 +142,7 @@ def from_root_file( # HACK: simple way of reading the number of active pixels within # the define ROI if roi is not None: + active_pixels = roi.active_pixels(nside_cnt, system="equatorial", ordering="RING") for pix_id in active_pixels: @@ -153,8 +158,14 @@ def from_root_file( data_tree_prefix = f"nHit{name}/data/count" bkg_tree_prefix = f"nHit{name}/bkg/count" - counts = map_infile[data_tree_prefix].array().to_numpy() - bkg = map_infile[bkg_tree_prefix].array().to_numpy() + # if using number of transitions different from the maptree + # then make sure counts are scaled accordingly. + counts: np.ndarray = map_infile[data_tree_prefix].array().to_numpy() * ( + n_transits / max(n_durations) + ) + bkg: np.ndarray = map_infile[bkg_tree_prefix].array().to_numpy() * ( + n_transits / max(n_durations) + ) except uproot.KeyInFileError: From 13d60521ad0b0f36127852c781cf4a3c29ced406 Mon Sep 17 00:00:00 2001 From: Ramiro Torres-Escobedo Date: Fri, 16 Sep 2022 10:43:19 -0500 Subject: [PATCH 3/6] scaling counts in obs. and bkg. for HDF maptrees --- hawc_hal/maptree/from_hdf5_file.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index 00bca4c..acbaa37 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/hawc_hal/maptree/from_hdf5_file.py @@ -1,6 +1,7 @@ from __future__ import absolute_import import collections +from curses import meta from hawc_hal.region_of_interest import get_roi_from_dict from hawc_hal.serialize import Serialization @@ -73,6 +74,8 @@ def from_hdf5_file(map_tree_file, roi, transits): data_analysis_bins = collections.OrderedDict() if transits is not None: + # if using user-specified transits, make sure to scale the counts + # in observation map and background map accordingly. n_transits = transits else: n_transits: float = meta_df["n_transits"].max() @@ -89,12 +92,14 @@ def from_hdf5_file(map_tree_file, roi, transits): # Read only the pixels that the user wants observation_hpx_map = SparseHealpix( - this_df.loc[active_pixels_user, "observation"].values, + this_df.loc[active_pixels_user, "observation"].values + * (n_transits / meta_df["n_transits"].max()), active_pixels_user, this_meta["nside"], ) background_hpx_map = SparseHealpix( - this_df.loc[active_pixels_user, "background"].values, + this_df.loc[active_pixels_user, "background"].values + * (n_transits / meta_df["n_transits"].max()), active_pixels_user, this_meta["nside"], ) @@ -102,8 +107,12 @@ def from_hdf5_file(map_tree_file, roi, transits): else: # Full sky - observation_hpx_map = DenseHealpix(this_df.loc[:, "observation"].values) - background_hpx_map = DenseHealpix(this_df.loc[:, "background"].values) + observation_hpx_map = DenseHealpix( + this_df.loc[:, "observation"].values * (n_transits / meta_df["n_transits"].max()) + ) + background_hpx_map = DenseHealpix( + this_df.loc[:, "background"].values * (n_transits / meta_df["n_transits"].max()) + ) # This signals the DataAnalysisBin that we are dealing with a full sky map active_pixels_user = None From 38768a9c8804cbcfa53cb61b1d31f8cd3ea0e7fd Mon Sep 17 00:00:00 2001 From: Ramiro Torres-Escobedo Date: Tue, 20 Sep 2022 13:26:05 -0500 Subject: [PATCH 4/6] minor code refactoring --- hawc_hal/maptree/from_hdf5_file.py | 20 +++++++++++--------- hawc_hal/maptree/from_root_file.py | 25 +++++++++++++++---------- 2 files changed, 26 insertions(+), 19 deletions(-) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index acbaa37..fa0d1fa 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/hawc_hal/maptree/from_hdf5_file.py @@ -3,9 +3,10 @@ import collections from curses import meta +from threeML.io.logging import setup_logger + from hawc_hal.region_of_interest import get_roi_from_dict from hawc_hal.serialize import Serialization -from threeML.io.logging import setup_logger log = setup_logger(__name__) log.propagate = False @@ -73,12 +74,11 @@ def from_hdf5_file(map_tree_file, roi, transits): data_analysis_bins = collections.OrderedDict() - if transits is not None: - # if using user-specified transits, make sure to scale the counts - # in observation map and background map accordingly. - n_transits = transits - else: - n_transits: float = meta_df["n_transits"].max() + assert ( + transits <= meta_df["transits"].max() + ), "Cannot use higher value than that of maptree" + + n_transits = meta_df["n_transits"].max() if transits is None else transits for bin_name in bin_names: @@ -108,10 +108,12 @@ def from_hdf5_file(map_tree_file, roi, transits): # Full sky observation_hpx_map = DenseHealpix( - this_df.loc[:, "observation"].values * (n_transits / meta_df["n_transits"].max()) + this_df.loc[:, "observation"].values + * (n_transits / meta_df["n_transits"].max()) ) background_hpx_map = DenseHealpix( - this_df.loc[:, "background"].values * (n_transits / meta_df["n_transits"].max()) + this_df.loc[:, "background"].values + * (n_transits / meta_df["n_transits"].max()) ) # This signals the DataAnalysisBin that we are dealing with a full sky map diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index b6b90b5..9491a32 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/hawc_hal/maptree/from_root_file.py @@ -55,7 +55,7 @@ def from_root_file( # Check that they exists and can be read if not file_existing_and_readable(map_tree_file): # pragma: no cover - # raise IOError("MapTree %s does not exist or is not readable" % map_tree_file) + raise IOError(f"MapTree {map_tree_file} does not exist or is not readable") # Make sure we have a proper ROI (or None) @@ -117,13 +117,12 @@ def from_root_file( maptree_durations = map_infile["BinInfo/totalDuration"].array() n_durations: np.ndarray = np.divide(maptree_durations, 24.0) - if transits is not None: - - n_transits = transits + assert transits <= max( + n_durations + ), "Cannot use a higher value than that of maptree." - else: - - n_transits: float = max(n_durations) + # use value of maptree unless otherwise specified by user + n_transits = max(n_durations) if transits is None else transits n_bins: int = data_bins_labels.shape[0] nside_cnt: int = hp.pixelfunc.npix2nside(npix_cnt) @@ -131,7 +130,9 @@ def from_root_file( # so far, a value of Nside of 1024 (perhaps will change) and a # RING HEALPix - assert nside_cnt == nside_bkg, "Nside value needs to be the same for counts and bkg. maps" + assert ( + nside_cnt == nside_bkg + ), "Nside value needs to be the same for counts and bkg. maps" assert scheme == 0, "NESTED HEALPix is not currently supported." @@ -143,7 +144,9 @@ def from_root_file( # the define ROI if roi is not None: - active_pixels = roi.active_pixels(nside_cnt, system="equatorial", ordering="RING") + active_pixels = roi.active_pixels( + nside_cnt, system="equatorial", ordering="RING" + ) for pix_id in active_pixels: @@ -185,7 +188,9 @@ def from_root_file( counts_hpx = SparseHealpix( counts[healpix_map_active > 0], active_pixels, nside_cnt ) - bkg_hpx = SparseHealpix(bkg[healpix_map_active > 0], active_pixels, nside_cnt) + bkg_hpx = SparseHealpix( + bkg[healpix_map_active > 0], active_pixels, nside_cnt + ) this_data_analysis_bin = DataAnalysisBin( name, From f7aaa0c6ee48bb879f6625dbc0adcbb5bdf2c4f3 Mon Sep 17 00:00:00 2001 From: Ramiro Torres-Escobedo Date: Tue, 20 Sep 2022 13:58:55 -0500 Subject: [PATCH 5/6] minor fix --- hawc_hal/maptree/from_hdf5_file.py | 6 +++--- hawc_hal/maptree/from_root_file.py | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index fa0d1fa..7087679 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/hawc_hal/maptree/from_hdf5_file.py @@ -74,12 +74,12 @@ def from_hdf5_file(map_tree_file, roi, transits): data_analysis_bins = collections.OrderedDict() + n_transits = meta_df["n_transits"].max() if transits is None else transits + assert ( - transits <= meta_df["transits"].max() + n_transits <= meta_df["transits"].max() ), "Cannot use higher value than that of maptree" - n_transits = meta_df["n_transits"].max() if transits is None else transits - for bin_name in bin_names: this_df = analysis_bins_df.loc[bin_name] diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index 9491a32..aed90fd 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/hawc_hal/maptree/from_root_file.py @@ -117,13 +117,13 @@ def from_root_file( maptree_durations = map_infile["BinInfo/totalDuration"].array() n_durations: np.ndarray = np.divide(maptree_durations, 24.0) - assert transits <= max( - n_durations - ), "Cannot use a higher value than that of maptree." - # use value of maptree unless otherwise specified by user n_transits = max(n_durations) if transits is None else transits + assert n_transits <= max( + n_durations + ), "Cannot use a higher value than that of maptree." + n_bins: int = data_bins_labels.shape[0] nside_cnt: int = hp.pixelfunc.npix2nside(npix_cnt) nside_bkg: int = hp.pixelfunc.npix2nside(npix_bkg) From be0f72163809a35e7d48042cd81d507c37a3aaf0 Mon Sep 17 00:00:00 2001 From: Ramiro Torres-Escobedo Date: Tue, 20 Sep 2022 14:12:04 -0500 Subject: [PATCH 6/6] minor bug fix --- hawc_hal/maptree/from_hdf5_file.py | 6 +++--- hawc_hal/maptree/from_root_file.py | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/hawc_hal/maptree/from_hdf5_file.py b/hawc_hal/maptree/from_hdf5_file.py index 7087679..0431621 100644 --- a/hawc_hal/maptree/from_hdf5_file.py +++ b/hawc_hal/maptree/from_hdf5_file.py @@ -76,9 +76,9 @@ def from_hdf5_file(map_tree_file, roi, transits): n_transits = meta_df["n_transits"].max() if transits is None else transits - assert ( - n_transits <= meta_df["transits"].max() - ), "Cannot use higher value than that of maptree" + # assert ( + # n_transits <= meta_df["transits"].max() + # ), "Cannot use higher value than that of maptree" for bin_name in bin_names: diff --git a/hawc_hal/maptree/from_root_file.py b/hawc_hal/maptree/from_root_file.py index aed90fd..c9b7622 100644 --- a/hawc_hal/maptree/from_root_file.py +++ b/hawc_hal/maptree/from_root_file.py @@ -120,9 +120,9 @@ def from_root_file( # use value of maptree unless otherwise specified by user n_transits = max(n_durations) if transits is None else transits - assert n_transits <= max( - n_durations - ), "Cannot use a higher value than that of maptree." + # assert n_transits <= max( + # n_durations + # ), "Cannot use a higher value than that of maptree." n_bins: int = data_bins_labels.shape[0] nside_cnt: int = hp.pixelfunc.npix2nside(npix_cnt)