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

@@ -35,25 +35,69 @@ 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,
+        ),
+        (
+            "big",
+            [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,
+        ),
+        (
+            "bottom",
+            [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,
+        ),
+        (
+            "bottom+small",
+            [0.01766961, 0.01693795, 0.01623237, 0.01987508] * u.ct / u.keV / u.s,
+            [0.16637264, 0.15948359, 0.15284002, 0.18713889] * u.ct / u.keV / u.s / u.cm**2,
+            [0.18145868, 0.18359402, 0.18161734, 0.17643195] * 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
     )
 
 

stixpy/calibration/visibility.py

@@ -94,6 +94,7 @@
     time_range: Time,
     energy_range: Quantity["energy"],  # noqa: F821
     flare_location: SkyCoord,
+    pixels: str = "big",
     no_shadowing: bool = False,
 ):
     r"""
@@ -109,6 +110,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 +170,26 @@
 
     e_cor_high, e_cor_low = get_elut_correction(e_ind, pixel_data)
 
+    # Get counts and other data.
+    pixel_slices = {
+        "all": slice(None),
+        "big": slice(0, 8),
+        "small": slice(8, None),
+        "top": slice(0, 4),
+        "bottom": slice(4, 8),
+        "bottom+small": slice(4, None),
+    }
+    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 +204,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 +220,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 +302,20 @@
     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["pixels"]
+    if pix_set in {"big", "top", "bottom"}:
+        phase += 46.1 * u.deg  # Center of large pixel in terms morie pattern phase
+    elif pix_set == "small":
+        phase += 22.5 * u.deg
+    elif pix_set == "all":
+        raise NotImplementedError("phase correction for combined big and small pixels not yet implemented.")
+    else:
+        raise ValueError("Set of pixels used to make meta pixels not recognised: {pix_set}")
 
     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))