Enable Meta Pixels to be generated from just top or bottom pixels

changelog/148.feature.rst

@@ -0,0 +1 @@
+Enable `stixpy.calibration.visibility.create_meta_pixels` to include only top or bottom row of big pixels.

stixpy/calibration/tests/test_visibility.py

@@ -3,8 +3,10 @@
 import pytest
 from astropy.tests.helper import assert_quantity_allclose
 from astropy.time import Time
+from numpy.ma.testutils import assert_equal
+from sunpy.time import TimeRange
 
-from stixpy.calibration.visibility import create_meta_pixels, get_uv_points_data
+from stixpy.calibration.visibility import create_meta_pixels, create_visibility, get_uv_points_data
 from stixpy.coordinates.frames import STIXImaging
 from stixpy.product import Product
 
@@ -35,25 +37,63 @@ def test_get_uv_points_data():
     assert uv_data["label"][0] == "3c"
 
 
-def test_create_meta_pixels(background_cpd):
+@pytest.mark.parametrize(
+    "pixel_set,expected_abcd_rate_kev,expected_abcd_rate_kev_cm0,expected_abcd_rate_kev_cm1",
+    [
+        (
+            "all",
+            [0.03509001, 0.03432438, 0.03248172, 0.03821136] * u.ct / u.keV / u.s,
+            [0.17336964, 0.16958685, 0.1604828, 0.1887913] * u.ct / u.keV / u.s / u.cm**2,
+            [0.18532299, 0.17855843, 0.1802053, 0.17609046] * u.ct / u.keV / u.s / u.cm**2,
+        ),
+        (
+            "top+bot",
+            [0.0339154, 0.03319087, 0.03131242, 0.03684958] * u.ct / u.keV / u.s,
+            [0.17628464, 0.17251869, 0.16275495, 0.19153585] * u.ct / u.keV / u.s / u.cm**2,
+            [0.18701911, 0.18205339, 0.18328145, 0.17945563] * u.ct / u.keV / u.s / u.cm**2,
+        ),
+        (
+            "small",
+            [0.00117461, 0.00113351, 0.00116929, 0.00136178] * u.ct / u.keV / u.s,
+            [0.1173439, 0.11323753, 0.11681252, 0.13604156] * u.ct / u.keV / u.s / u.cm**2,
+            [0.15272384, 0.11138607, 0.12108232, 0.11141279] * u.ct / u.keV / u.s / u.cm**2,
+        ),
+        (
+            "top",
+            [0.01742041, 0.01738642, 0.01624934, 0.01833627] * u.ct / u.keV / u.s,
+            [0.18109474, 0.18074145, 0.16892087, 0.19061566] * u.ct / u.keV / u.s / u.cm**2,
+            [0.18958941, 0.17299885, 0.17864632, 0.17571344] * u.ct / u.keV / u.s / u.cm**2,
+        ),
+        (
+            "bot",
+            [0.01649499, 0.01580445, 0.01506308, 0.01851331] * u.ct / u.keV / u.s,
+            [0.17147454, 0.16429592, 0.15658903, 0.19245605] * u.ct / u.keV / u.s / u.cm**2,
+            [0.18444881, 0.19110794, 0.18791658, 0.18319781] * u.ct / u.keV / u.s / u.cm**2,
+        ),
+    ],
+)
+def test_create_meta_pixels(
+    background_cpd, pixel_set, expected_abcd_rate_kev, expected_abcd_rate_kev_cm0, expected_abcd_rate_kev_cm1
+):
     time_range = Time(["2022-08-24T14:00:37.271", "2022-08-24T14:50:17.271"])
     energy_range = [20, 76] * u.keV
     meta_pixels = create_meta_pixels(
         background_cpd,
         time_range=time_range,
         energy_range=energy_range,
         flare_location=STIXImaging(0 * u.arcsec, 0 * u.arcsec),
+        pixels=pixel_set,
         no_shadowing=True,
     )
 
+    assert_quantity_allclose(expected_abcd_rate_kev, meta_pixels["abcd_rate_kev"][0, :], atol=1e-7 * u.ct / u.keV / u.s)
+
     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, :],
+        expected_abcd_rate_kev_cm0, meta_pixels["abcd_rate_kev_cm"][0, :], atol=1e-7 * expected_abcd_rate_kev_cm0.unit
     )
 
     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, :],
+        expected_abcd_rate_kev_cm1, meta_pixels["abcd_rate_kev_cm"][-1, :], atol=1e-7 * expected_abcd_rate_kev_cm1.unit
     )
 
 
@@ -95,3 +135,29 @@ def test_create_meta_pixels_timebins(flare_cpd):
     )
 
     assert_quantity_allclose(np.sum(flare_cpd.duration[0:3]), meta_pixels["time_range"].dt.to(u.s))
+
+
+@pytest.mark.parametrize(
+    "pix_set, real_comp",
+    [
+        ("blahblah", None),
+        ("top+bot", np.cos(46.1 * u.deg)),
+        ("all", np.cos(45 * u.deg)),
+        ("small", np.cos(22.5 * u.deg)),
+    ],
+)
+def test_create_visibility(pix_set, real_comp):
+    # counts chosen to make real component 0 so real comp will only be phase from pixel combination
+    fake_meta_pixels = {
+        "abcd_rate_kev_cm": np.repeat([[0, 0, 1, 0]], 32, axis=0) * u.ct,
+        "abcd_rate_error_kev_cm": np.repeat([[0, 0, 0, 0]], 32, axis=0) * u.ct,
+        "energy_range": [4, 10] * u.keV,
+        "time_range": TimeRange("2025-01-30", "2025-01-31"),
+        "pixels": pix_set,
+    }
+    if pix_set == "blahblah":
+        with pytest.raises(ValueError):
+            create_visibility(fake_meta_pixels)
+    else:
+        vis = create_visibility(fake_meta_pixels)
+        assert_equal(np.real(vis[0].visibilities.value), real_comp)

stixpy/calibration/visibility.py

@@ -10,13 +10,15 @@
 from astropy.units import Quantity
 from sunpy.coordinates import HeliographicStonyhurst
 from sunpy.time import TimeRange
+from xrayvision.visibility import Visibilities, VisMeta
 
 from stixpy.calibration.energy import get_elut
 from stixpy.calibration.grid import get_grid_transmission
 from stixpy.calibration.livetime import get_livetime_fraction
 from stixpy.coordinates.frames import STIXImaging
 from stixpy.coordinates.transforms import get_hpc_info
 from stixpy.io.readers import read_subc_params
+from stixpy.product.sources import CompressedPixelData, RawPixelData, SummedCompressedPixelData
 
 __all__ = [
     "get_subcollimator_info",
@@ -26,9 +28,13 @@
     "calibrate_visibility",
 ]
 
-from xrayvision.visibility import Visibilities, VisMeta
-
-from stixpy.product.sources import CompressedPixelData, RawPixelData, SummedCompressedPixelData
+_PIXEL_SLICES = {
+    "top": slice(0, 4),
+    "bot": slice(4, 8),
+    "top+bot": slice(0, 8),
+    "all": slice(None),
+    "small": slice(8, None),
+}
 
 
 def get_subcollimator_info():
@@ -94,6 +100,7 @@
     time_range: Time,
     energy_range: Quantity["energy"],  # noqa: F821
     flare_location: SkyCoord,
+    pixels: str = "top+bot",
     no_shadowing: bool = False,
 ):
     r"""
@@ -109,6 +116,9 @@
         Start and end energies
     flare_location
         The coordinates of flare used to calculate grid transmission
+    pixels :
+        The set of pixels to use to create the meta pixels.
+        Allowed values are 'all', 'big' (default), 'small', 'top', 'bottom'.
     no_shadowing : bool optional
         If set to True turn grid shadowing correction off
     Returns
@@ -166,14 +176,18 @@
 
     e_cor_high, e_cor_low = get_elut_correction(e_ind, pixel_data)
 
+    # Get counts and other data.
+    idx_pix = _PIXEL_SLICES.get(pixels.lower(), None)
+    if idx_pix is None:
+        raise ValueError(f"Unrecognised input for 'pixels': {pixels}. Supported values: {list(_PIXEL_SLICES.keys())}")
     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]
+    ct = counts[t_ind][..., idx_pix, e_ind]
+    ct[..., 0] = ct[..., 0] * e_cor_low[..., idx_pix]
+    ct[..., -1] = ct[..., -1] * e_cor_high[..., idx_pix]
+    ct_error = count_errors[t_ind][..., idx_pix, e_ind]
+    ct_error[..., 0] = ct_error[..., 0] * e_cor_low[..., idx_pix]
+    ct_error[..., -1] = ct_error[..., -1] * e_cor_high[..., idx_pix]
 
     lt = (livefrac * pixel_data.data["timedel"].reshape(-1, 1).to("s"))[t_ind].sum(axis=0)
 
@@ -188,10 +202,8 @@
         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_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_counts = ct_summed.reshape(ct_summed.shape[0], -1, 4).sum(axis=1)
+    abcd_count_errors = np.sqrt((ct_error_summed.reshape(ct_error_summed.shape[0], -1, 4) ** 2).sum(axis=1))
 
     abcd_rate = abcd_counts / lt.reshape(-1, 1)
     abcd_rate_error = abcd_count_errors / lt.reshape(-1, 1)
@@ -206,13 +218,15 @@
                    0.096194997, 0.096194997, 0.010009999, 0.010009999, 0.010009999, 0.010009999,] * u.cm**2
     # fmt: on
 
-    areas = pixel_areas.reshape(-1, 4)[0:2].sum(axis=0)
+    areas = pixel_areas[idx_pix].reshape(-1, 4).sum(axis=0)
 
     meta_pixels = {
+        "abcd_rate_kev": abcd_rate_kev,
         "abcd_rate_kev_cm": abcd_rate_kev / areas,
         "abcd_rate_error_kev_cm": abcd_rate_error_kev / areas,
         "time_range": time_range,
         "energy_range": energy_range,
+        "pixels": pixels,
     }
 
     return meta_pixels
@@ -286,21 +300,22 @@
     amplitude_error = np.sqrt((real / observed_amplitude * real_err) ** 2 + (imag / observed_amplitude * imag_err) ** 2)
 
     # Apply pixel correction
-    phase += 46.1 * u.deg  # Center of large pixel in terms morie pattern phase
-    # TODO add correction for small pixel
+    pix_set = meta_pixels.get("pixels", None)
+    if pix_set in {"top", "bot", "top+bot"}:
+        phase += 46.1 * u.deg  # Center of large pixel in terms morie pattern phase
+    elif pix_set == "all":
+        phase += 45 * u.deg
+    elif pix_set == "small":
+        phase += 22.5 * u.deg
+    else:
+        raise ValueError(
+            f"Set of pixels used to make meta pixels not recognised, {pix_set} should be one of {list(_PIXEL_SLICES.keys())}"
+        )
 
     obsvis = (np.cos(phase) + np.sin(phase) * 1j) * observed_amplitude
 
     imaging_detector_indices = vis_info["isc"] - 1
 
-    # vis = SimpleNamespace(
-    #     obsvis=obsvis[imaging_detector_indices],
-    #     amplitude=observed_amplitude[imaging_detector_indices],
-    #     amplitude_error=amplitude_error[imaging_detector_indices],
-    #     phase=phase[imaging_detector_indices],
-    #     **vis_info,
-    # )
-
     vis_meta = VisMeta(
         instrumet="STIX",
         spectral_range=meta_pixels["energy_range"],

stixpy/coordinates/tests/test_transforms.py

@@ -62,7 +62,7 @@ def test_stx_to_hpc_obstime_end():
     stix_coord_rt = hp.transform_to(STIXImaging(obstime=times[0], obstime_end=times[-1]))
 
     stix_coord_rt_interp = hp.transform_to(STIXImaging(obstime=times[1]))  # noqa: F841
-    assert_quantity_allclose(0 * u.deg, stix_coord.separation(stix_coord_rt), atol=1e-17 * u.deg)
+    assert_quantity_allclose(0 * u.deg, stix_coord.separation(stix_coord_rt), atol=1e-9 * u.deg)
     assert np.all(stix_coord.obstime.isclose(stix_coord_rt.obstime))
     assert np.all(stix_coord.obstime_end.isclose(stix_coord_rt.obstime_end))