Clean up and add test

changelog/106.bugfix.rst

@@ -0,0 +1 @@
+Fix ELUT correction when the requested data has less then than 32 energy channels.

changelog/106.feature.rst

@@ -1 +1 @@
-Implement background subtraction in the demo, fix ELUT correction, implement ELUT correction for background counts, implement Energy Edge Masks
+Add background subtraction step to imaging demo, update `EnergyEdgeMasks`.

examples/imaging_demo.py

@@ -53,8 +53,7 @@
 # Set time and energy ranges which will be considered for the science and the background file
 
 time_range_sci = ["2021-09-23T15:21:00", "2021-09-23T15:24:00"]
-# time_range_bkg = ["2021-09-23T09:59:23.757", "2021-09-23T11:35:23.757"] # Set this range larger than the actual observation time
-time_range_bkg = ["2021-09-23T09:00:00", "2021-09-23T12:00:00"]
+time_range_bkg = ["2021-09-23T09:00:00", "2021-09-23T12:00:00"]  # Set this range larger than the actual observation time
 energy_range = [28, 40]
 
 ###############################################################################

stixpy/calibration/tests/test_visibility.py

@@ -1,11 +1,41 @@
 import astropy.units as u
+import pytest
+from astropy.tests.helper import assert_quantity_allclose
 
-from stixpy.calibration.visibility import get_uv_points_data
+from stixpy.calibration.visibility import create_meta_pixels, get_uv_points_data
+from stixpy.product import Product
+
+
+@pytest.fixture
+@pytest.mark.remote_data
+def flare_cpd():
+    return Product('https://pub099.cs.technik.fhnw.ch/fits/L1/2022/08/28/SCI/solo_L1_stix-sci-xray-cpd_20220828T154401-20220828T161600_V02_2208284257-61808.fits')
+
+
+@pytest.fixture
+@pytest.mark.remote_data
+def background_cpd():
+    return Product('https://pub099.cs.technik.fhnw.ch/fits/L1/2022/08/24/SCI/solo_L1_stix-sci-xray-cpd_20220824T140037-20220824T145017_V02_2208242079-61784.fits')
 
 
 def test_get_uv_points_data():
     uv_data = get_uv_points_data()
-    assert uv_data['u'][0]     == -0.03333102271709666 / u.arcsec
-    assert uv_data['v'][0]     == 0.005908224739704219 / u.arcsec
-    assert uv_data['isc'][0]   == 1
+    assert uv_data['u'][0] == -0.03333102271709666 / u.arcsec
+    assert uv_data['v'][0] == 0.005908224739704219 / u.arcsec
+    assert uv_data['isc'][0] == 1
     assert uv_data['label'][0] == '3c'
+
+
+def test_create_meta_pixels(background_cpd):
+    time_range = ['2022-08-24T14:00:37.271', '2022-08-24T14:50:17.271']
+    energy_range = [20, 76]
+    meta_pixels = create_meta_pixels(background_cpd, time_range=time_range, energy_range=energy_range,
+                                     phase_center=[0,0]*u.arcsec, no_shadowing=True)
+
+    assert_quantity_allclose([0.17628464, 0.17251869,
+                              0.16275495, 0.19153585] * u.ct / (u.keV * u.cm ** 2 * u.s),
+                             meta_pixels['abcd_rate_kev_cm'][0, :])
+
+    assert_quantity_allclose([0.18701911, 0.18205339,
+                              0.18328145, 0.17945563] * u.ct / (u.keV * u.cm ** 2 * u.s),
+                             meta_pixels['abcd_rate_kev_cm'][-1, :])

stixpy/calibration/visibility.py

@@ -88,7 +88,7 @@ def create_meta_pixels(pixel_data, time_range, energy_range, phase_center, no_sh
     energy_range
         Start and end energies
     phase_center
-        The coordinates of the phase center
+        The coordinates of the center
     no_shadowing : bool optional
         If set to True turn grid shadowing correction off
     Returns
@@ -157,32 +157,7 @@ def create_meta_pixels(pixel_data, time_range, energy_range, phase_center, no_sh
 
     pixel_data.data["livefrac"] = livefrac
 
-    energy_mask = pixel_data.energy_masks.energy_mask.astype(bool)
-
-    elut = get_elut(pixel_data.time_range.center)
-    ebin_edges_low = np.zeros((32,12,32), dtype=float)
-    ebin_edges_low[...,1:] = elut.e_actual
-    ebin_edges_low = ebin_edges_low[...,energy_mask]
-    ebin_edges_high = np.zeros((32, 12, 32), dtype=float)
-    ebin_edges_high[...,0:-1] = elut.e_actual
-    ebin_edges_high[...,-1] = np.nan
-    ebin_edges_high = ebin_edges_high[..., energy_mask]
-
-    ebin_widths = ebin_edges_high - ebin_edges_low
-
-    ebin_sci_edges_low = elut.e[..., 0:-1].value
-    ebin_sci_edges_low = ebin_sci_edges_low[..., energy_mask]
-    ebin_sci_edges_high = elut.e[..., 1:].value
-    ebin_sci_edges_high = ebin_sci_edges_high[..., energy_mask]
-
-    e_cor_low = (ebin_edges_high[..., e_ind[0]] - ebin_sci_edges_low[..., e_ind[0]]) / ebin_widths[..., e_ind[0]]
-    e_cor_high = (ebin_sci_edges_high[..., e_ind[-1]] - ebin_edges_low[..., e_ind[-1]]) / ebin_widths[..., e_ind[-1]]
-
-    e_cor = (  # noqa
-        (ebin_edges_high[..., e_ind[-1]] - ebin_edges_low[..., e_ind[0]])
-        * u.keV
-        / (elut.e[e_ind[-1] + 1] - elut.e[e_ind[0]])
-    )
+    e_cor_high, e_cor_low = get_elut_correction(e_ind, pixel_data)
 
     counts = pixel_data.data["counts"].astype(float)
     count_errors = np.sqrt(pixel_data.data["counts_comp_err"].astype(float).value**2 + counts.value) * u.ct
@@ -236,6 +211,42 @@ def create_meta_pixels(pixel_data, time_range, energy_range, phase_center, no_sh
     return meta_pixels
 
 
+def get_elut_correction(e_ind, pixel_data):
+    r"""
+    Get ELUT correction factors
+
+    Only correct the low and high energy edges as internal bins are contiguous.
+
+    Parameters
+    ----------
+    e_ind : `np.ndarray`
+        Energy channel indices
+    pixel_data : `~stixpy.products.sources.CompressedPixelData`
+        Pixel data
+
+    Returns
+    -------
+
+    """
+    energy_mask = pixel_data.energy_masks.energy_mask.astype(bool)
+    elut = get_elut(pixel_data.time_range.center)
+    ebin_edges_low = np.zeros((32, 12, 32), dtype=float)
+    ebin_edges_low[..., 1:] = elut.e_actual
+    ebin_edges_low = ebin_edges_low[..., energy_mask]
+    ebin_edges_high = np.zeros((32, 12, 32), dtype=float)
+    ebin_edges_high[..., 0:-1] = elut.e_actual
+    ebin_edges_high[..., -1] = np.nan
+    ebin_edges_high = ebin_edges_high[..., energy_mask]
+    ebin_widths = ebin_edges_high - ebin_edges_low
+    ebin_sci_edges_low = elut.e[..., 0:-1].value
+    ebin_sci_edges_low = ebin_sci_edges_low[..., energy_mask]
+    ebin_sci_edges_high = elut.e[..., 1:].value
+    ebin_sci_edges_high = ebin_sci_edges_high[..., energy_mask]
+    e_cor_low = (ebin_edges_high[..., e_ind[0]] - ebin_sci_edges_low[..., e_ind[0]]) / ebin_widths[..., e_ind[0]]
+    e_cor_high = (ebin_sci_edges_high[..., e_ind[-1]] - ebin_edges_low[..., e_ind[-1]]) / ebin_widths[..., e_ind[-1]]
+    return e_cor_high, e_cor_low
+
+
 def create_visibility(meta_pixels):
     r"""
     Create visibilities from meta-pixels
@@ -286,9 +297,9 @@ def create_visibility(meta_pixels):
 
 
 @u.quantity_input
-def get_uv_points_data(d_det: u.Quantity[u.mm]=47.78 * u.mm, d_sep:u.Quantity[u.mm]=545.30 * u.mm):
+def get_uv_points_data(d_det: u.Quantity[u.mm] = 47.78 * u.mm, d_sep:u.Quantity[u.mm] = 545.30 * u.mm):
     r"""
-    Returns the STIX (u,v) points coordinates defined in [1]. The coordinates are ordered with respect to the detector index.
+    Return the STIX (u,v) points coordinates defined in [1], ordered with respect to the detector index.
 
     Parameters
     ----------

stixpy/map/stix.py

@@ -17,14 +17,14 @@ def plot(self, **kwargs):
 
         # Until this is resolved need to manually override https://github.com/astropy/astropy/issues/16339
         res = super().plot(**kwargs)
+        if self.coordinate_system == 'stiximaging':
+            res.axes.coords[0].set_coord_type('longitude', coord_wrap=180*u.deg)
+            res.axes.coords[1].set_coord_type('latitude')
 
-        res.axes.coords[0].set_coord_type('longitude', coord_wrap=180*u.deg)
-        res.axes.coords[1].set_coord_type('latitude')
+            res.axes.coords[0].set_format_unit(u.arcsec)
+            res.axes.coords[1].set_format_unit(u.arcsec)
 
-        res.axes.coords[0].set_format_unit(u.arcsec)
-        res.axes.coords[1].set_format_unit(u.arcsec)
-
-        res.axes.coords[0].set_axislabel(r'STIX $\Theta_{X}$')
-        res.axes.coords[1].set_axislabel(r'STIX $\Theta_{Y}$')
+            res.axes.coords[0].set_axislabel(r'STIX $\Theta_{X}$')
+            res.axes.coords[1].set_axislabel(r'STIX $\Theta_{Y}$')
 
         return res