Feat calibration fixes (#106)

Calibration fixes - Fixed ELUT correction - Implemented ELUT correction for count erros - Defined Energy Edge Masks - Implemented background subtraction in the demo --------- Co-authored-by: Shane Maloney <shane.maloney@dias.ie>
TCDSolar · May 3, 2024 · a35a42c · a35a42c
1 parent d44a1c8
commit a35a42c
Show file tree

Hide file tree

Showing 8 changed files with 162 additions and 67 deletions.
diff --git a/changelog/106.bugfix.rst b/changelog/106.bugfix.rst
@@ -0,0 +1 @@
+Fix ELUT correction when the requested data has less then than 32 energy channels.
diff --git a/changelog/106.feature.rst b/changelog/106.feature.rst
@@ -0,0 +1 @@
+Add background subtraction step to imaging demo, update `EnergyEdgeMasks`.
diff --git a/docs/tutorials/quickstart.rst b/docs/tutorials/quickstart.rst
@@ -171,7 +171,7 @@ will return the correct product type. In this case we are providing the path to
     PixelMasks
     [0...339]: [['1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1' '1']]
 
-    EnergyMasks
+    EnergyEdgeMasks
     [0]: [0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]
 
 A spectrogram plot can be obtained by call the `plot_spectrogram` method.
@@ -226,7 +226,7 @@ Now let look at the pixel data, the data can be loaded just same as the spectrog
     PixelMasks
     [0...146]: [['1' '1' '1' '1' '1' '1' '1' '1' '0' '0' '0' '0']]
 
-    EnergyMasks
+    EnergyEdgeMasks
     [0]: [_,_,2,3,4,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_,_]
 
 Pixel data supports the same plot methods as spectrogram but notice the plots only contain

diff --git a/examples/imaging_demo.py b/examples/imaging_demo.py
@@ -34,35 +34,55 @@
 logger.setLevel("DEBUG")
 
 #############################################################################
-# Create a product
+# Read science file as Product
 
-cpd = Product(
+cpd_sci = Product(
     "http://pub099.cs.technik.fhnw.ch/fits/L1/2021/09/23/SCI/solo_L1_stix-sci-xray-cpd_20210923T152015-20210923T152639_V02_2109230030-62447.fits"
 )
-cpd
+cpd_sci
+
+#############################################################################
+# Read background file as Product
+
+cpd_bkg = Product(
+    "http://pub099.cs.technik.fhnw.ch/fits/L1/2021/09/23/SCI/solo_L1_stix-sci-xray-cpd_20210923T095923-20210923T113523_V02_2109230083-57078.fits"
+)
+cpd_bkg
 
 ###############################################################################
-# Set time and energy ranges which will be used to create the visibilties
+# Set time and energy ranges which will be considered for the science and the background file
 
-time_range = ["2021-09-23T15:21:00", "2021-09-23T15:24:00"]
+time_range_sci = ["2021-09-23T15:21:00", "2021-09-23T15:24: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]
 
 ###############################################################################
-# Creat the meta pixel, A, B, C, D
+# Create the meta pixel, A, B, C, D for the science and the background data
+
+meta_pixels_sci = create_meta_pixels(
+    cpd_sci, time_range=time_range_sci, energy_range=energy_range, phase_center=[0, 0] * u.arcsec, no_shadowing=True
+)
 
-meta_pixels = create_meta_pixels(
-    cpd, time_range=time_range, energy_range=energy_range, phase_center=[0, 0] * u.arcsec, no_shadowing=True
+meta_pixels_bkg = create_meta_pixels(
+    cpd_bkg, time_range=time_range_bkg, energy_range=energy_range, phase_center=[0, 0] * u.arcsec, no_shadowing=True
 )
 
+###############################################################################
+# Perform background subtraction
+
+meta_pixels_bkg_subtracted = {"abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
+                              "abcd_rate_error_kev_cm": np.sqrt( meta_pixels_sci["abcd_rate_error_kev_cm"]**2 +
+                                                                 meta_pixels_bkg["abcd_rate_error_kev_cm"]**2)}
+
 ###############################################################################
 # Create visibilites from the meta pixels
 
-vis = create_visibility(meta_pixels)
+vis = create_visibility(meta_pixels_bkg_subtracted)
 
 ###############################################################################
-# Calibrate the visibilties
+# Calibrate the visibilities
 
-# Extra phase calihraiton not needed with these
+# Extra phase calibration not needed with these
 uv_data = get_uv_points_data()
 vis.u = uv_data['u']
 vis.v = uv_data['v']
@@ -95,7 +115,7 @@
 # Make a full disk back projection (inverse transform) map
 
 fd_bp_map = vis_to_map(stix_vis, imsize, pixel_size=pixel)
-fd_bp_map.meta["rsun_obs"] = cpd.meta["rsun_arc"]
+fd_bp_map.meta["rsun_obs"] = cpd_sci.meta["rsun_arc"]
 fd_bp_map.draw_limb()
 fd_bp_map.plot()
 
@@ -113,16 +133,20 @@
 axes.plot_coord(max_hpc, marker=".", markersize=50, fillstyle="none", color="r", markeredgewidth=2)
 axes.set_xlabel("STIX Y")
 axes.set_ylabel("STIX X")
-axes.set_title(f'STIX {" ".join(time_range)} {"-".join([str(e) for e in energy_range])} keV')
+axes.set_title(f'STIX {" ".join(time_range_sci)} {"-".join([str(e) for e in energy_range])} keV')
 
 ################################################################################
 # Use estimated flare location to create more accurate visibilities
 
-meta_pixels = create_meta_pixels(
-    cpd, time_range=time_range, energy_range=energy_range, phase_center=[max_hpc.Tx, max_hpc.Ty], no_shadowing=False
+meta_pixels_sci = create_meta_pixels(
+    cpd_sci, time_range=time_range_sci, energy_range=energy_range, phase_center=[max_hpc.Tx, max_hpc.Ty], no_shadowing=True
 )
 
-vis = create_visibility(meta_pixels)
+meta_pixels_bkg_subtracted = {"abcd_rate_kev_cm": meta_pixels_sci["abcd_rate_kev_cm"] - meta_pixels_bkg["abcd_rate_kev_cm"],
+                              "abcd_rate_error_kev_cm": np.sqrt( meta_pixels_sci["abcd_rate_error_kev_cm"]**2 +
+                                                                 meta_pixels_bkg["abcd_rate_error_kev_cm"]**2)}
+
+vis = create_visibility(meta_pixels_bkg_subtracted)
 uv_data = get_uv_points_data()
 vis.u = uv_data['u']
 vis.v = uv_data['v']
@@ -202,7 +226,7 @@
         setattr(stix_vis1, k, v[idx])
 
 em_map = em(
-    meta_pixels["abcd_rate_kev_cm"],
+    meta_pixels_bkg_subtracted["abcd_rate_kev_cm"],
     stix_vis1,
     shape=imsize,
     pixel_size=pixel,

diff --git a/stixpy/calibration/tests/test_visibility.py b/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, :])
diff --git a/stixpy/calibration/visibility.py b/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,38 +157,29 @@ def create_meta_pixels(pixel_data, time_range, energy_range, phase_center, no_sh
 
     pixel_data.data["livefrac"] = livefrac
 
-    elut = get_elut(pixel_data.time_range.center)
-    ebin_edges_low = elut.e_actual[..., 0:-1]
-    ebin_edges_high = elut.e_actual[..., 1:]
-    ebin_widths = ebin_edges_high - ebin_edges_low
-
-    e_cor_low = (ebin_edges_high[..., e_ind[0]] - elut.e[..., e_ind[0] + 1].value) / ebin_widths[..., e_ind[0]]
-    e_cor_high = (elut.e[..., e_ind[-1] + 2].value - ebin_edges_high[..., e_ind[-1] - 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"]
-    ct = counts[t_ind][..., e_ind].astype(float)
+    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
+    ct = counts[t_ind][..., e_ind]
     ct[..., 0:8, 0] = ct[..., 0:8, 0] * e_cor_low[..., 0:8]
     ct[..., 0:8, -1] = ct[..., 0:8, -1] * e_cor_high[..., 0:8]
+    ct_error = count_errors[t_ind][..., e_ind]
+    ct_error[..., 0:8, 0] = ct_error[..., 0:8, 0] * e_cor_low[..., 0:8]
+    ct_error[..., 0:8, -1] = ct_error[..., 0:8, -1] * e_cor_high[..., 0:8]
 
     lt = (livefrac * pixel_data.data["timedel"].reshape(-1, 1).to("s"))[t_ind].sum(axis=0)
 
-    ct_sumed = ct.sum(axis=(0, 3)).astype(float)
-
-    err_sumed = (np.sqrt(ct.sum(axis=(0, 3)).value)) * u.ct
+    ct_summed = ct.sum(axis=(0, 3))#.astype(float)
+    ct_error_summed = np.sqrt(np.sum(ct_error**2, axis=(0, 3)))
 
     if not no_shadowing:
         grid_shadowing = get_grid_transmission(phase_center)
-        ct_sumed = ct_sumed / grid_shadowing.reshape(-1, 1) / 4  # transmission grid ~ 0.5*0.5 = .25
-        err_sumed = err_sumed / grid_shadowing.reshape(-1, 1) / 4
+        ct_summed = ct_summed / grid_shadowing.reshape(-1, 1) / 4  # transmission grid ~ 0.5*0.5 = .25
+        ct_error_summed = ct_error_summed / grid_shadowing.reshape(-1, 1) / 4
 
-    abcd_counts = ct_sumed.reshape(ct_sumed.shape[0], -1, 4)[:, [0, 1], :].sum(axis=1)
-    abcd_count_errors = np.sqrt((err_sumed.reshape(ct_sumed.shape[0], -1, 4)[:, [0, 1], :] ** 2).sum(axis=1))
+    abcd_counts = ct_summed.reshape(ct_summed.shape[0], -1, 4)[:, [0, 1], :].sum(axis=1)
+    abcd_count_errors = np.sqrt((ct_error_summed.reshape(ct_error_summed.shape[0], -1, 4)[:, [0, 1], :] ** 2).sum(axis=1))
 
     abcd_rate = abcd_counts / lt.reshape(-1, 1)
     abcd_rate_error = abcd_count_errors / lt.reshape(-1, 1)
@@ -213,13 +204,49 @@ def create_meta_pixels(pixel_data, time_range, energy_range, phase_center, no_sh
         0.010009999,
     ] * u.cm**2
 
-    areas = np.array(pixel_areas).reshape(-1, 4)[0:2].sum(axis=0)
+    areas = pixel_areas.reshape(-1, 4)[0:2].sum(axis=0)
 
     meta_pixels = {"abcd_rate_kev_cm": abcd_rate_kev / areas, "abcd_rate_error_kev_cm": abcd_rate_error_kev / areas}
 
     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
@@ -270,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
     ----------

diff --git a/stixpy/map/stix.py b/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