Merge pull request #339 from israelmcmc/PR307_merge_v0.3

PR307 and v0.3.x merge conflict resolution. Follow up on "Initial patches to make cosipy work with histpy 2.x"

Author: israelmcmc

cosipy/response/FullDetectorResponse.py

@@ -9,7 +9,7 @@
 import copy
 
 class SourceInjector():
-    
+
     def __init__(self, response_path, response_frame = "spacecraftframe"):
 
         """
@@ -24,23 +24,23 @@ def __init__(self, response_path, response_frame = "spacecraftframe"):
         """
 
         self.response_path =  response_path
-        
+
         if response_frame == "spacecraftframe" or response_frame == "galactic":
-            
+
             self.response_frame = response_frame
-        
+
         else:
             raise ValueError("The response frame can only be `spacecraftframe` or `galactic`!")
 
 
 
     @staticmethod
     def get_psr_in_galactic(coordinate, response_path, spectrum):
-        
+
         """
         Get the point source response (psr) in galactic. Please be aware that you must use a galactic response!
         To do: to make the weight parameter not hardcoded
-        
+
         Parameters
         ----------
         coordinate : astropy.coordinates.SkyCoord
@@ -49,39 +49,27 @@ def get_psr_in_galactic(coordinate, response_path, spectrum):
             The path to the response.
         spectrum : astromodels.functions
             The spectrum of the source to be placed at the hypothesis coordinate.
-        
+
         Returns
         -------
         psr : histpy.Histogram
             The point source response of the spectrum at the hypothesis coordinate.
         """
-        
-        # Open the response
-        # Notes from Israel: Inside it contains a single histogram with all the regular axes for a Compton Data Space (CDS) analysis, in galactic coordinates. Since there is no class yet to handle it, this is how to read in the HDF5 manually.
-        
-        with h5.File(response_path) as f:
 
-            axes_group = f['hist/axes']
-            axes = []
-            for axis in axes_group.values():
-                # Get class. Backwards compatible with version
-                # with only Axis
-                axis_cls = Axis
-                if '__class__' in axis.attrs:
-                    class_module, class_name = axis.attrs['__class__']
-                    axis_cls = getattr(sys.modules[class_module], class_name)
-                axes += [axis_cls._open(axis)]
-        axes = Axes(axes)
-        
-        # get the pixel number of the hypothesis coordinate
-        map_temp = HealpixMap(base = axes[0])
-        coordinate_pix_number = map_temp.ang2pix(coordinate)
-        
-        # get the expectation for the hypothesis coordinate (a point source)
+        # Open the response Histogram, but don't read the whole thing into memory --
+        # just get the axes, the specific contents pixel(s), and its unit
         with h5.File(response_path) as f:
+
+            axes = Axes.open(f['hist/axes'])
+
+            # get the pixel number of the hypothesis coordinate
+            hp_axis = axes[0]
+            coordinate_pix_number = hp_axis.ang2pix(coordinate)
+
+            # get the expectation for the hypothesis coordinate (a point source)
             pix = coordinate_pix_number
-            psr = PointSourceResponse(axes[1:], f['hist/contents'][pix+1], unit = f['hist'].attrs['unit'])
-                
+            psr = PointSourceResponse(axes[1:], f['hist/contents'][pix], unit = f['hist'].attrs['unit'])
+
         return psr
 
 
@@ -107,26 +95,28 @@ def inject_point_source(self, spectrum, coordinate, orientation = None, source_n
             The path to save the injected data to a `.h5` file. This should include the file name. (the default is `None`, which means the injected data won't be saved.
         project_axes : list, optional
             The axes to project before saving the data file (the default is `None`, which means the data won't be projected).
-            
+
         Returns
         -------
         histpy.Histogram
             The `Histogram object of the injected spectrum.`
         """
-        
-        
+
+
         # get the point source response in local frame
         if self.response_frame == "spacecraftframe":
 
             if orientation == None:
                 raise TypeError("The when the data are binned in spacecraftframe frame, orientation must be provided to compute the expected counts.")
 
+
+
             with FullDetectorResponse.open(self.response_path) as response:
-                
+
                 scatt_map = orientation.get_scatt_map(coordinate, response.nside*2, coordsys = 'galactic', earth_occ = True)
-                
+
                 psr = response.get_point_source_response(coord=coordinate, scatt_map=scatt_map)
-                
+
         # get the point source response in galactic frame
         elif self.response_frame == "galactic":
 
@@ -223,62 +213,61 @@ def inject_extended_source(
             injected.write(data_save_path)
 
         return injected
-    
+
     def inject_model(self, model, orientation = None, make_spectrum_plot = False, data_save_path = None, project_axes = None):
-        
+
         if self.response_frame == "spacecraftframe":
             if orientation == None:
                 raise TypeError("The when the data are binned in spacecraftframe frame, orientation must be provided to compute the expected counts.")
-                
+
         self.components = {}
-        
+
         # first inject point sources
         point_sources = model.point_sources
-        
+
         # iterate through all point sources
         for name, source in point_sources.items():
-            
+
             injected = self.inject_point_source(spectrum = source.spectrum.main.shape, coordinate = source.position.sky_coord,
                                                 orientation = orientation, source_name = name)
-            
+
             injected.axes["Em"].axis_scale = "log"  # set to log scale manually. This inconsistency is from the detector response module
-            
+
             self.components[name] = injected
-            
+
         # second inject extended sources
         extended_sources = model.extended_sources
-        
+
         # iterate through all extended sources
         for name, source in extended_sources.items():
-            
+
             injected = self.inject_extended_source(source_model = source, source_name = name)
             self.components[name] = injected
-            
-        
+
+
         if len(self.components) == 1:
-            
+
             # if there is only one component, the injected all is just the only component
             injected_all = list(self.components.values())[0]
-            
+
             if data_save_path is not None:
                 injected_all.write(data_save_path)
-                
+
         elif len(self.components) > 1:
-            
+
             injected_list = list(self.components.values())
-            
+
             injected_all = copy.deepcopy(injected_list[0])
-            
+
             # add the rest of the injected sources
             for i in injected_list[1:]:
                 injected_all += i
-                
+
         if make_spectrum_plot:
             ax, plot = injected_all.project("Em").draw(color="green", linewidth=2)
             ax.set_xscale("log")
             ax.set_yscale("log")
             ax.set_xlabel("Em [keV]", fontsize=14, fontweight="bold")
             ax.set_ylabel("Counts", fontsize=14, fontweight="bold")
-                
+
             return injected_all
-        

cosipy/source_injector/source_injector.py

@@ -4,5 +4,5 @@ SIMPLE  =                    T / conforms to FITS standard
 P
 �
 �
-
@@k�O�z�@1ѝ�t@D�U��-K@R8K�#�?������?������?�s@?������?������?�s@?������?������?�s@?������@k����^�@k���Q�@k�b�x�5@k�(�:�@k��U��,@k���)�<@k�x1m��@k�=ʙ��@k����b@k��T2�@0�1�m$@0�`�wv@0�Ta��Z@1
- 4Wq~@1�ح�@1%�Nh�@@13��Q8@1AM�,"*@1O��]�@1\�N��@D�c�{f@D�w��F@D����4�@D��KD�)@D��V��@D��Ԧ<�@E
�Ŷw@E�*fW�@E?��@E%P�X8@RG���@RD]�΢$@R@�痓j@R=w��#�@R:	Խ�@R6�,e�0@R3Z���@R/����u@R,9ڏ�@R(�,J�
+
@@k�O�z�@1ѝ�t@D�U��-L@R8K�#�?������?������?�s@?������?������?�s@?������?������?�s@?������@k����^�@k���Q�@k�b�x�5@k�(�:�@k��U��,@k���)�<@k�x1m��@k�=ʙ��@k����b@k��T2�@0�1�m%@0�`�wv@0�Ta��Z@1
+ 4Wq~@1�ح�@1%�Nh�@@13��Q8@1AM�,"*@1O��]�@1\�N��@D�c�{d@D�w��F@D����4�@D��KD�)@D��V��@D��Ԧ<�@E
�Ŷw@E�*fW�@E?��@E%P�X8@RG���@RD]�΢#@R@�痓j@R=w��#�@R:	Խ�@R6�,e�0@R3Z���@R/����u@R,9ڏ�@R(�,J�