diff --git a/changelog/148.feature.rst b/changelog/148.feature.rst new file mode 100644 index 0000000..0b4e367 --- /dev/null +++ b/changelog/148.feature.rst @@ -0,0 +1 @@ +Enable `stixpy.calibration.visibility.create_meta_pixels` to include only top or bottom row of big pixels. diff --git a/stixpy/calibration/tests/test_visibility.py b/stixpy/calibration/tests/test_visibility.py index 9550f5e..606d757 100644 --- a/stixpy/calibration/tests/test_visibility.py +++ b/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,7 +37,44 @@ 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( @@ -43,17 +82,18 @@ def test_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) diff --git a/stixpy/calibration/visibility.py b/stixpy/calibration/visibility.py index f88dedf..c990ca2 100644 --- a/stixpy/calibration/visibility.py +++ b/stixpy/calibration/visibility.py @@ -10,6 +10,7 @@ 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 @@ -17,6 +18,7 @@ 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 @@ def create_meta_pixels( 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 @@ def create_meta_pixels( 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 @@ def create_meta_pixels( 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 @@ def create_meta_pixels( 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 @@ def create_meta_pixels( 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 @@ def create_visibility(meta_pixels): 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"], diff --git a/stixpy/coordinates/tests/test_transforms.py b/stixpy/coordinates/tests/test_transforms.py index 7c0072a..f66ef65 100644 --- a/stixpy/coordinates/tests/test_transforms.py +++ b/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))