diff --git a/scopesim/detector/detector.py b/scopesim/detector/detector.py index b99c5dba..9e133621 100644 --- a/scopesim/detector/detector.py +++ b/scopesim/detector/detector.py @@ -25,6 +25,7 @@ class Detector: def __init__(self, header, cmds=None, **kwargs): + self._native_header = header.copy() self._hdu = ImageHDU(header=header, data=zeros_from_header(header)) self.meta = {} self.meta.update(header) @@ -51,10 +52,9 @@ def extract_from(self, image_plane, spline_order=1, reset=True): self.header.update(sky_wcs.to_header()) def reset(self): - """Reset internal HDU data to all-zeros-array.""" - # The detector might have been converted to integers by the - # Quantization effect, so it is not possible to use zeros_like. - self._hdu.data = np.zeros(self._hdu.data.shape) + """Reset charge and restore the configured, unbinned detector geometry.""" + self._hdu = ImageHDU(header=self._native_header, data=zeros_from_header(self._native_header)) + self.meta.update(self._native_header) @property def hdu(self): diff --git a/scopesim/effects/detector_list.py b/scopesim/effects/detector_list.py index 7cb07514..4a3c5ea5 100644 --- a/scopesim/effects/detector_list.py +++ b/scopesim/effects/detector_list.py @@ -330,16 +330,17 @@ def detector_headers(self, ids=None): hdr["PC2_1"] = -sang hdr["PC2_2"] = cang - # hdr["GAIN"] = row["gain"] - if "id" in row: - hdr["DET_ID"] = row["id"] - hdr["EXTNAME"] = f"DET_{row['id']}" - row_dict = { f"HIERARCH {col.upper()}" if len(col) > 8 else col: row[col] for col in row.colnames } hdr.update(row_dict) + if "id" in row.colnames: + # Keep routing identifiers numeric when detector metadata is + # serialized to FITS (NumPy integers otherwise become strings). + hdr["ID"] = int(row["id"]) + hdr["DET_ID"] = int(row["id"]) + hdr["EXTNAME"] = f"DET_{row['id']}" hdrs.append(hdr) return hdrs diff --git a/scopesim/effects/electronic/__init__.py b/scopesim/effects/electronic/__init__.py index c35dfe06..e086896c 100644 --- a/scopesim/effects/electronic/__init__.py +++ b/scopesim/effects/electronic/__init__.py @@ -52,9 +52,9 @@ def apply_to(self, obj, **kwargs): from .electrons import LinearityCurve, ADConversion, InterPixelCapacitance from .noise import (Bias, PoorMansHxRGReadoutNoise, BasicReadoutNoise, - PixelResponseNonUniformity, ShotNoise, DarkCurrent) + PixelResponseNonUniformity, ShotNoise, DarkCurrent, ReadoutNoiseFromTable) from .exposure import AutoExposure, ExposureIntegration, ExposureOutput -from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage +from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage, OnChipBinning from .dmps import DetectorModePropertiesSetter diff --git a/scopesim/effects/electronic/dmps.py b/scopesim/effects/electronic/dmps.py index 296703a8..2cf3e213 100644 --- a/scopesim/effects/electronic/dmps.py +++ b/scopesim/effects/electronic/dmps.py @@ -23,6 +23,9 @@ class DetectorModePropertiesSetter(Effect): mode_properties : dict A dictionary containing the DET parameters to be changed for each mode. See below for an example yaml entry. + mode_key : str + Command selecting the mode, by default ``!OBS.detector_readout_mode``. + Separate keys allow independent detector arms to select their modes. Examples -------- @@ -70,21 +73,20 @@ def __init__(self, **kwargs): check_keys(self.meta, self.required_keys, action="error") self.mode_properties = kwargs["mode_properties"] + self.mode_key = kwargs.get("mode_key", "!OBS.detector_readout_mode") def apply_to(self, obj, **kwargs): logger.debug("Entering DetectorModePropertiesSetter.apply_to") - mode_name = kwargs.get( - "detector_readout_mode", - from_currsys("!OBS.detector_readout_mode", self.cmds), - ) - + mode_name = kwargs.get("detector_readout_mode", + from_currsys(self.mode_key, + self.cmds)) if isinstance(obj, ImagePlane) and mode_name == "auto": mode_name = self.select_mode(obj, **kwargs) logger.info("Detector mode set to %s", mode_name) self.meta["detector_readout_mode"] = mode_name props_dict = self.mode_properties[mode_name] - self.cmds["!OBS.detector_readout_mode"] = mode_name + self.cmds[self.mode_key] = mode_name for key, value in props_dict.items(): self.cmds[key] = value logger.debug("DMPS setting %s: %s", key, value) diff --git a/scopesim/effects/electronic/electrons.py b/scopesim/effects/electronic/electrons.py index c2bee25b..ac3e1323 100644 --- a/scopesim/effects/electronic/electrons.py +++ b/scopesim/effects/electronic/electrons.py @@ -263,7 +263,7 @@ def __init__(self, **kwargs): super().__init__(**kwargs) params = { "dtype": "uint16", - "gain": 1. # default, usually overridden from yaml + "gain": "!DET.gain" } self.meta.update(params) self.meta.update(kwargs) @@ -297,17 +297,18 @@ def _get_gain(self, det_id: int) -> Real: # Apply the gain value (copy from DarkCurrent) # Note that this does not cater for the case where the gain is given # as a plain dictionary. Should we implement that? - if hasattr(self.cmds["!DET.gain"], "dic"): - gain = self.cmds["!DET.gain"].dic[det_id] + configured_gain = from_currsys(self.meta["gain"], self.cmds) + if hasattr(configured_gain, "dic"): + gain = configured_gain.dic[det_id] logger.info(f"Detector {det_id}: applying gain {gain}") return gain - if isinstance(self.cmds["!DET.gain"], Real): - gain = self.cmds["!DET.gain"] + if isinstance(configured_gain, Real): + gain = configured_gain logger.info(f"Applying gain {gain}") return gain raise ValueError( f"{self.__class__.__name__}.meta['gain'] must be either " - f"dict or float, but is {self.cmds['!DET.gain']}") + f"dict or float, but is {configured_gain}") def _apply_to_det(self, det: Detector) -> None: """Subclasses can override if more params needed in call.""" diff --git a/scopesim/effects/electronic/exposure.py b/scopesim/effects/electronic/exposure.py index 342e9ac6..f1123737 100644 --- a/scopesim/effects/electronic/exposure.py +++ b/scopesim/effects/electronic/exposure.py @@ -287,6 +287,9 @@ def apply_to(self, obj, **kwargs): "readout call." ) + obj.header["EXPTIME"] = (dit * ndit, "Total accumulated integration time [s]") + obj.header['DIT'] = (dit, 'Integration time per exposure [s]') + obj.header['NDIT'] = ndit obj.data = obj.data * dit * ndit return obj diff --git a/scopesim/effects/electronic/noise.py b/scopesim/effects/electronic/noise.py index 55521854..8c6aef4f 100644 --- a/scopesim/effects/electronic/noise.py +++ b/scopesim/effects/electronic/noise.py @@ -54,7 +54,7 @@ def __call__( return data + dark_level * dit * ndit def _get_dark_level(self, det_id: int) -> float: - dark_level = float(from_currsys(self.meta["value"], self.cmds)) + dark_level = from_currsys(self.meta["value"], self.cmds) if isinstance(dark_level, Real): return dark_level if isinstance(dark_level, Mapping): @@ -286,6 +286,23 @@ def _pseudo_random_field( return image +class ReadoutNoiseFromTable(BasicReadoutNoise): + """Interpolate measured read noise in electrons as a function of DIT [s].""" + + required_keys = {"filename", "dit", "ndit"} + + def apply_to(self, det, **kwargs): + if not isinstance(det, Detector): + return det + dit = from_currsys(self.meta["dit"], self.cmds) + times = self.table["dit"] + if not times[0] <= dit <= times[-1]: + raise ValueError(f"DIT {dit} s is outside the measured read-noise table") + self.meta["noise_std"] = np.interp(dit, times, self.table["read_noise"]) + det.header["RDNOISE"] = (self.meta["noise_std"], "Read noise per integration [electron]") + return super().apply_to(det, **kwargs) + + class PixelResponseNonUniformity(RandomEffect): """Pixel Response Non-Uniformity (PRNU). diff --git a/scopesim/effects/electronic/pixels.py b/scopesim/effects/electronic/pixels.py index e26b26e9..efb966cf 100644 --- a/scopesim/effects/electronic/pixels.py +++ b/scopesim/effects/electronic/pixels.py @@ -93,19 +93,55 @@ class BinnedImageBase(ElectronicEffect): z_order: ClassVar[tuple[int, ...]] = (870,) + def __call__(self, data: ArrayLike) -> NDArray: + bx, by = self.bin_factors + ny, nx = data.shape + return data.reshape(ny//by, by, nx//bx, bx).sum(axis=(1, 3)) + + def bin_header(self, native_header): + """WCS of the summed pixel blocks, also used by trace coordinates.""" + bx, by = self.bin_factors + if not isinstance(bx, int) or not isinstance(by, int) or min(bx, by) < 1: + raise ValueError('Bin factors must be positive integers') + if 'allowed_bins' in self.meta and [bx, by] not in self.meta['allowed_bins']: + raise ValueError(f'Unsupported binning: {bx} x {by}') + header = native_header.copy() + for axis, factor in ((1, bx), (2, by)): + if header[f'NAXIS{axis}'] % factor: + raise ValueError('Binning must divide the full detector dimensions') + header[f'NAXIS{axis}'] //= factor + for suffix in [''] + list('ABCDEFGHIJKLMNOPQRSTUVWXYZ'): + if f'CRPIX1{suffix}' not in header: + continue + for axis, factor in ((1, bx), (2, by)): + header[f'CRPIX{axis}{suffix}'] = (header[f'CRPIX{axis}{suffix}'] - .5)/factor + .5 + if f'CD1_1{suffix}' in header: + for i in (1, 2): + for j, factor in ((1, bx), (2, by)): + header[f'CD{i}_{j}{suffix}'] *= factor + else: + header[f'CDELT1{suffix}'] *= bx + header[f'CDELT2{suffix}'] *= by + if f'PC1_2{suffix}' in header: + header[f'PC1_2{suffix}'] *= by/bx + if f'PC2_1{suffix}' in header: + header[f'PC2_1{suffix}'] *= bx/by + header['BINX'], header['BINY'] = bx, by + return header + + def _apply_to_det(self, det): + header = self.bin_header(det.header) + det.data = self(det.data) + det.header.update(header) + det.meta.update(header) + class BinnedImage(BinnedImageBase): required_keys = {"bin_size"} - def __call__(self, data: ArrayLike) -> NDArray: - height, width = data.shape - binned_data = data.reshape(( - height//self.bin_size, - self.bin_size, - width//self.bin_size, - self.bin_size - )).sum(axis=3).sum(axis=1) - return binned_data + @property + def bin_factors(self): + return self.bin_size, self.bin_size @property def bin_size(self) -> int: @@ -115,15 +151,9 @@ def bin_size(self) -> int: class UnequalBinnedImage(BinnedImageBase): required_keys = {"binx","biny"} - def __call__(self, data: ArrayLike) -> NDArray: - height, width = data.shape - binned_data = data.reshape(( - height//self.biny, - self.biny, - width//self.binx, - self.binx - )).sum(axis=3).sum(axis=1) - return binned_data + @property + def bin_factors(self): + return self.binx, self.biny @property def binx(self) -> int: @@ -132,3 +162,14 @@ def binx(self) -> int: @property def biny(self) -> int: return from_currsys(self.meta["biny"], self.cmds) + + +class OnChipBinning(UnequalBinnedImage): + """Sum charge in FITS x/y pixel blocks before adding read noise. + + Unlike a software sum after readout, each output pixel receives one read + noise draw. Configure this effect between charge integration and read noise. + WCS reference pixels refer to the centers of the new pixel blocks. + """ + + z_order = (815,) diff --git a/scopesim/effects/fits_headers.py b/scopesim/effects/fits_headers.py index 9110f333..987117ed 100644 --- a/scopesim/effects/fits_headers.py +++ b/scopesim/effects/fits_headers.py @@ -153,6 +153,9 @@ class ExtraFitsKeywords(Effect): header['HIERARCH ESO DET DIT'] = (5, '[s] exposure length') header['EXTNAME'] = "DET2.DATA" + A group may additionally specify ``detector_id``. It then applies only to + HDULists containing that ``DET_ID``, including their primary headers. + Resolved and un-resolved keywords --------------------------------- ScopeSim uses bang-strings to resolve global parameters. @@ -285,6 +288,12 @@ def apply_to(self, hdul, **kwargs): return hdul for dic in self.dict_list: + # Each detector array produces its own HDUList. A detector selector + # also permits arm-specific metadata in that array's primary HDU. + if "detector_id" in dic: + detector_ids = [hdu.header["DET_ID"] for hdu in hdul if "DET_ID" in hdu.header] + if dic["detector_id"] not in detector_ids: + continue resolved = flatten_dict(dic.get("keywords", {}), resolve=True, optics_manager=opt_train) unresolved = flatten_dict(dic.get("unresolved_keywords", {})) diff --git a/scopesim/effects/mosaic_trace_list.py b/scopesim/effects/mosaic_trace_list.py index 34e6bb66..c74c9d21 100644 --- a/scopesim/effects/mosaic_trace_list.py +++ b/scopesim/effects/mosaic_trace_list.py @@ -255,6 +255,9 @@ def apply_to(self, det, **kwargs): image = det.data detwcs = WCS(det.header, key="D") + # ExposureIntegration has already recorded the actual accumulated time. + # New table HDUs need those cards, but not the detector-image WCS. + exposure_header = fits.Header(card for card in det.header.cards if card.keyword in ("DIT", "NDIT", "EXPTIME")) output_format = from_currsys(self.meta["format"], self.cmds) @@ -309,6 +312,7 @@ def apply_to(self, det, **kwargs): else: logger.error("Unknown output format: %s", output_format) + det.header.update(exposure_header) return det diff --git a/scopesim/optics/optical_train.py b/scopesim/optics/optical_train.py index cd8ab65b..6974bbe9 100644 --- a/scopesim/optics/optical_train.py +++ b/scopesim/optics/optical_train.py @@ -122,6 +122,7 @@ def __init__(self, cmds=None): self.detector_managers = [] self.yaml_dicts = None self._last_source = None + self._readout_id = 0 if cmds is not None: self.load(cmds) @@ -234,12 +235,13 @@ def trace_detector_coordinates(self, *, xi=0 * u.arcsec, explicit definition of which output image owns each trace coordinate. """ from ..effects import ( - BinnedImage, + DetectorModePropertiesSetter, + ImagePlaneBackground, Rotate90CCD, SelectorWheel, SpectralTraceList, - UnequalBinnedImage, ) + from ..effects.electronic.pixels import BinnedImageBase trace_lists = self.optics_manager.get_all(SpectralTraceList) if len(trace_lists) != 1: @@ -248,14 +250,18 @@ def trace_detector_coordinates(self, *, xi=0 * u.arcsec, "SpectralTraceList.") trace_list = trace_lists[0] - if self.optics_manager.image_plane_effects: - raise NotImplementedError( - "Trace detector coordinates do not yet support image-plane " - "effects.") - if self.optics_manager.detector_array_effects: - raise NotImplementedError( - "Trace detector coordinates do not yet support detector-array " - "effects.") + for effect in self.optics_manager.image_plane_effects: + selected = effect.wheel_effects.values() if isinstance(effect, SelectorWheel) else [effect] + if any(not isinstance(item, ImagePlaneBackground) for item in selected): + raise NotImplementedError('Trace detector coordinates only support additive image-plane backgrounds.') + for effect in self.optics_manager.detector_array_effects: + if not isinstance(effect, DetectorModePropertiesSetter): + raise NotImplementedError('Trace detector coordinates do not yet support this detector-array effect.') + # Mode setters update commands, not image coordinates. Resolve an + # explicit mode before reading any binning command it can set. + if from_currsys(effect.mode_key, self.cmds) == 'auto': + raise NotImplementedError('Select an explicit detector readout mode before inspecting coordinates.') + effect.apply_to(None) fovs_by_trace = {} for fov in self.fov_manager.fovs: @@ -326,31 +332,14 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): "detectors on one image plane.") detector = detector_manager[0] detector_id = detector.header["ID"] + readout_header = detector._native_header.copy() for effect in self.optics_manager.detector_effects: - selected_effect = ( - effect.wheel_effects.get(detector_id) - if isinstance(effect, SelectorWheel) - else effect - ) + selected_effect = effect.get_effect(detector_id) if isinstance(effect, SelectorWheel) else effect if selected_effect is None: continue - if isinstance(selected_effect, BinnedImage): - bin_size = from_currsys( - selected_effect.meta["bin_size"], self.cmds) - if bin_size != 1: - raise NotImplementedError( - "Trace detector coordinates do not yet support " - "post-extraction detector binning.") - if isinstance(selected_effect, UnequalBinnedImage): - binx = from_currsys( - selected_effect.meta["binx"], self.cmds) - biny = from_currsys( - selected_effect.meta["biny"], self.cmds) - if binx != 1 or biny != 1: - raise NotImplementedError( - "Trace detector coordinates do not yet support " - "post-extraction detector binning.") + if isinstance(selected_effect, BinnedImageBase): + readout_header = selected_effect.bin_header(readout_header) if isinstance(selected_effect, Rotate90CCD): rotations = from_currsys( selected_effect.meta["rotations"], self.cmds) @@ -359,7 +348,7 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): "Trace detector coordinates do not yet support " "post-extraction detector rotation.") - detector_wcs = WCS(detector.header, key="D") + detector_wcs = WCS(detector._native_header, key="D") image_plane_wcs = image_plane._det_wcs if ( not np.allclose(image_plane_wcs.wcs.pc, np.eye(2)) @@ -429,13 +418,7 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): f"Trace {trace_name!r} has overlapping FOV wavelength " "ranges.") - ( - _, - _, - x_grid_mm, - y_grid_mm, - subimage_wcs, - ) = trace._focal_plane_grid(fov) + _, _, x_grid_mm, y_grid_mm, subimage_wcs = trace._focal_plane_grid(fov) if ( not np.allclose(subimage_wcs.wcs.pc, np.eye(2)) or not np.allclose( @@ -450,16 +433,7 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): x_mm = trace.xilam2x(xi_flat[in_fov], wave_flat[in_fov]) y_mm = trace.xilam2y(xi_flat[in_fov], wave_flat[in_fov]) - subimage_x = ( - (x_mm - x_grid_mm[0]) - / (x_grid_mm[-1] - x_grid_mm[0]) - * (len(x_grid_mm) - 1) - ) - subimage_y = ( - (y_mm - y_grid_mm[0]) - / (y_grid_mm[-1] - y_grid_mm[0]) - * (len(y_grid_mm) - 1) - ) + subimage_x, subimage_y = subimage_wcs.all_world2pix(x_mm, y_mm, 0) subimage_shape = np.array( [len(x_grid_mm), len(y_grid_mm)]) @@ -477,6 +451,8 @@ def image_origin_on_canvas(shape_xy, image_wcs, canvas_wcs): f"Requested wavelengths for trace {trace_name!r} " "do not map through a configured FOV.") + focal_x, focal_y = detector_wcs.all_pix2world(x_pix, y_pix, 0) + x_pix, y_pix = WCS(readout_header, key='D').all_world2pix(focal_x, focal_y, 0) point_count = wave_grid.size columns["readout_index"].extend( np.full(point_count, readout_index)) @@ -534,6 +510,7 @@ def observe(self, orig_source=None, update=True, **kwargs): """ if update: self.update(**kwargs) + self._readout_id = 0 if self.cmds["!SIM.random.seed"] is None: # Concretize seed from system entropy to keep constant during the @@ -731,7 +708,7 @@ def readout( self, filename: Path | str | None = None, reset: bool = True, - roid: int = 0, + roid: int | None = None, **kwargs ) -> list[HDUList]: """ @@ -746,13 +723,10 @@ def readout( reset : bool, optional If true (the default), do not keep any parameters set during any previous readout. This is usually what you want. - roid : int, optional - Readout ID, used for seeding randomness. If you run multiple - readouts on the same observation and want the random effects to - produce different results, set this parameter with e.g. the loop - counter. Simple integers are fine (0, 1, 2, ...), because this gets - wrangled into a proper random seed downstream. If you need to - reproduce a specific readout in a sequence, set e.g. `roid=4`. + roid : int or None, optional + Readout ID for reproducible noise. By default, successive calls + use 0, 1, 2, ... within each observation. Supply an explicit ID to + replay that noise realization with the same root seed. **kwargs : Any other kwargs passed to the detector effects. @@ -764,6 +738,9 @@ def readout( """ if reset: self.cmds.clear() + if roid is None: + roid = self._readout_id + self._readout_id += 1 # Hack to make sure AutoExposure and Quantization work properly. # Should probably be removed once #428 is fixed properly. @@ -774,7 +751,7 @@ def readout( # self.cmds["!OBS.ndit"] = None # TODO: This is still hacky but seems to work for now... params = {"exptime": None, "roid": roid} - params.update(kwargs) + params.update({key.removeprefix('!OBS.'): value for key, value in kwargs.items()}) if params["exptime"] is not None and params.get("dit") is None: params.update( dit=None, @@ -789,8 +766,8 @@ def readout( for key, value in params.items(): if is_bangkey(key): - # skip if already explicitly !ABC.xyz logger.debug("%s: %s given in kwargs", key, value) + self.cmds[key] = value continue # otherwise fallback to !OBS logger.debug("%s: %s given in kwargs, put in !OBS", key, value) diff --git a/scopesim/tests/tests_effects/test_OnChipBinning.py b/scopesim/tests/tests_effects/test_OnChipBinning.py new file mode 100644 index 00000000..ec3347f7 --- /dev/null +++ b/scopesim/tests/tests_effects/test_OnChipBinning.py @@ -0,0 +1,75 @@ +import numpy as np +import pytest +from scopesim.detector import Detector +from scopesim.optics.image_plane_utils import header_from_list_of_xy +from astropy.wcs import WCS +from scopesim.effects.electronic import OnChipBinning, UnequalBinnedImage, BasicReadoutNoise + + +def detector(nx=120, ny=80): + header = header_from_list_of_xy([-.015*nx/2, .015*nx/2], [-.015*ny/2, .015*ny/2], .015, 'D') + return Detector(header) + + + +@pytest.mark.parametrize('bx,by', [(1, 1), (1, 2), (2, 2)]) +def test_binning_charge_and_focal_coordinates(bx, by): + det = detector() + det.data[:] = 7 + original = WCS(det.header, key='D') + OnChipBinning(binx=bx, biny=by).apply_to(det) + assert det.data.shape == (80//by, 120//bx) + assert np.all(det.data == 7*bx*by) + assert det.data.sum() == 7*120*80 + binned = WCS(det.hdu.header, key='D') + expected = original.all_pix2world((np.array([10, 20])+.5)*bx-.5, (np.array([12, 23])+.5)*by-.5, 0) + np.testing.assert_allclose(binned.all_pix2world([10, 20], [12, 23], 0), expected, atol=1e-12) + det.reset() + assert det.data.shape == (80, 120) + assert det.header['CDELT2D'] == .015 + OnChipBinning(binx=bx, biny=by).apply_to(det) + assert det.data.shape == (80//by, 120//bx) + + + +def test_binning_rotated_cd_matrix(): + det = detector() + det.header['CD1_1D'], det.header['CD1_2D'] = .012, -.009 + det.header['CD2_1D'], det.header['CD2_2D'] = .009, .012 + wcs = WCS(det.header, key='D') + OnChipBinning(binx=1, biny=2).apply_to(det) + actual = WCS(det.header, key='D').all_pix2world(9, 15, 0) + np.testing.assert_allclose(actual, wcs.all_pix2world(9, 30.5, 0), atol=1e-12) + + + +def test_binning_rejects_unsupported_modes(): + with pytest.raises(ValueError, match='Unsupported'): + OnChipBinning(binx=2, biny=1, allowed_bins=[[1, 1], [1, 2], [2, 2]]).apply_to(detector()) + with pytest.raises(ValueError): + OnChipBinning(binx=1, biny=3).apply_to(detector()) + + + +def test_ccd_bin_does_not_multiply_amplifier_noise(): + det = detector(600, 400) + OnChipBinning(binx=2, biny=2).apply_to(det) + BasicReadoutNoise(noise_std=3, ndit=1, random_seed=17).apply_to(det) + assert det.data.std() == pytest.approx(3, rel=.02) + + +def test_software_sum_accumulates_independent_read_noise(): + det = detector(600, 400) + BasicReadoutNoise(noise_std=3, ndit=1, random_seed=17).apply_to(det) + UnequalBinnedImage(binx=2, biny=2).apply_to(det) + assert det.data.std() == pytest.approx(6, rel=.02) + assert det.header['CDELT1D'] == pytest.approx(.03) + assert det.header['CDELT2D'] == pytest.approx(.03) + + +def test_header_projection_matches_applied_binning(): + det = detector() + binning = OnChipBinning(binx=1, biny=2) + expected = binning.bin_header(det.header) + binning.apply_to(det) + assert det.header == expected diff --git a/scopesim/tests/tests_effects/test_ReadoutNoiseFromTable.py b/scopesim/tests/tests_effects/test_ReadoutNoiseFromTable.py new file mode 100644 index 00000000..ef372b48 --- /dev/null +++ b/scopesim/tests/tests_effects/test_ReadoutNoiseFromTable.py @@ -0,0 +1,25 @@ +import numpy as np +import pytest +from scopesim.detector import Detector +from scopesim.optics.image_plane_utils import header_from_list_of_xy +from scopesim.effects.electronic import ReadoutNoiseFromTable, ExposureOutput + + +def detector(nx=120, ny=80): + header = header_from_list_of_xy([-.015*nx/2, .015*nx/2], [-.015*ny/2, .015*ny/2], .015, 'D') + return Detector(header) + + + +@pytest.mark.parametrize('ndit', [1, 4, 9]) +def test_nir_averaging_preserves_signal_and_reduces_noise(tmp_path, ndit): + noise_file = tmp_path / 'noise.dat' + noise_file.write_text('dit read_noise\n10 12\n100 6\n') + det = detector(500, 400) + det.data[:] = 30*ndit + ReadoutNoiseFromTable(filename=str(noise_file), dit=55, ndit=ndit, random_seed=21).apply_to(det) + ExposureOutput(dit=55, ndit=ndit, mode='average').apply_to(det) + assert det.data.mean() == pytest.approx(30, abs=.1) + assert det.data.std() == pytest.approx(9/np.sqrt(ndit), rel=.02) + with pytest.raises(ValueError, match='outside'): + ReadoutNoiseFromTable(filename=str(noise_file), dit=101, ndit=ndit).apply_to(det) diff --git a/scopesim/tests/tests_effects/test_detector_arm_configuration.py b/scopesim/tests/tests_effects/test_detector_arm_configuration.py new file mode 100644 index 00000000..2b81134c --- /dev/null +++ b/scopesim/tests/tests_effects/test_detector_arm_configuration.py @@ -0,0 +1,53 @@ +import numpy as np +import pytest +from scopesim.detector import Detector +from scopesim.optics.image_plane_utils import header_from_list_of_xy +from astropy.io import fits +from scopesim.commands import UserCommands +from scopesim.effects import DetectorList, ExtraFitsKeywords, DetectorModePropertiesSetter +from scopesim.effects.electronic import ADConversion + + +def detector(nx=120, ny=80): + header = header_from_list_of_xy([-.015*nx/2, .015*nx/2], [-.015*ny/2, .015*ny/2], .015, 'D') + return Detector(header) + + + +def test_gain_uses_selected_arm_and_missing_gain_fails(): + det = detector() + det.data[:] = 12 + ADConversion(gain='!DET.gain_uvb', dtype='float32', cmds={'!DET.gain_uvb': .6}).apply_to(det) + np.testing.assert_allclose(det.data, 20) + with pytest.raises(ValueError, match='gain_uvb'): + ADConversion(gain='!DET.gain_uvb', dtype='float32', cmds={}).apply_to(det) + + + +def test_detector_specific_primary_headers(): + eff = ExtraFitsKeywords(header_dict=[{'detector_id': 1, 'ext_number': [0, 1], 'keywords': {'ARM': 'VIS'}}]) + for detector_id in [0, 1, 2]: + hdul = fits.HDUList([fits.PrimaryHDU(), fits.ImageHDU(header=fits.Header({'DET_ID': detector_id}))]) + eff.apply_to(hdul) + assert ('ARM' in hdul[0].header) == (detector_id == 1) + + + +def test_detector_layout_sets_real_id_header(): + effect = DetectorList(array_dict={'id': [2], 'x_cen': [0.], 'y_cen': [0.], 'x_size': [16], 'y_size': [8], 'pixel_size': [.018], 'angle': [0.], 'gain': [2.29]}, x_cen_unit='mm', y_cen_unit='mm', x_size_unit='pixel', y_size_unit='pixel', pixel_size_unit='mm', angle_unit='deg', gain_unit='electron/adu') + header = effect.detector_headers()[0] + assert header['DET_ID'] == 2 + assert header['EXTNAME'] == 'DET_2' + det = Detector(header) + assert det.meta['ID'] == 2 + assert det.hdu.header['DET_ID'] == 2 + + + +def test_independent_detector_readout_modes(): + cmds = UserCommands(properties={'!DET.mode_blue': 'slow', '!DET.mode_red': 'fast'}) + for arm in ['blue', 'red']: + modes = {'slow': {f'!DET.gain_{arm}': .6}, 'fast': {f'!DET.gain_{arm}': 1.8}} + DetectorModePropertiesSetter(mode_properties=modes, mode_key=f'!DET.mode_{arm}', cmds=cmds).apply_to(None) + assert cmds['!DET.gain_blue'] == .6 + assert cmds['!DET.gain_red'] == 1.8 diff --git a/scopesim/tests/tests_effects/test_mosaic_trace_list.py b/scopesim/tests/tests_effects/test_mosaic_trace_list.py index 39433af4..f5bef8e8 100644 --- a/scopesim/tests/tests_effects/test_mosaic_trace_list.py +++ b/scopesim/tests/tests_effects/test_mosaic_trace_list.py @@ -12,6 +12,7 @@ MosaicSpectralTraceList, MosaicSpectralTrace, MosaicOutputFormat) +from scopesim.effects.electronic.exposure import ExposureIntegration from scopesim.optics.image_plane_utils import header_from_list_of_xy from scopesim.detector import Detector @@ -201,3 +202,25 @@ def test_image_returns_detector(self, apply_image): def test_image_detector_has_image(self, apply_image): assert isinstance(apply_image._hdu, fits.ImageHDU) + + +@pytest.mark.usefixtures("patch_mock_path_mosaic") +@pytest.mark.parametrize("output_format", ["collapse1d", "table", "image"]) +def test_output_format_preserves_accumulated_counts_and_exposure(output_format): + header = header_from_list_of_xy([-25, 25], [-2048, 2048], 1, "D") + detector = Detector(header) + rows = [2023, 2031, 2039, 2047, 2055, 2063, 2071] + detector.data[rows, :] = 1 + rate_sum = detector.data.sum() + ExposureIntegration(dit=17., ndit=3).apply_to(detector) + output = MosaicOutputFormat(filename="TRACE_MOS-LR-J.fits", format=output_format).apply_to(detector) + assert output.header["DIT"] == 17. + assert output.header["NDIT"] == 3 + assert output.header["EXPTIME"] == 51. + if output_format == "image": + counts = output.data.sum() + assert "CTYPE1D" in output.header + else: + counts = output.data["spectrum"].sum() + assert "CTYPE1D" not in output.header + assert counts == pytest.approx(rate_sum * 51.) diff --git a/scopesim/tests/tests_optics/test_OpticalTrain.py b/scopesim/tests/tests_optics/test_OpticalTrain.py index c1931499..1585e350 100644 --- a/scopesim/tests/tests_optics/test_OpticalTrain.py +++ b/scopesim/tests/tests_optics/test_OpticalTrain.py @@ -23,6 +23,11 @@ Effect, SpectralTraceList, UnequalBinnedImage, + OnChipBinning, + BasicReadoutNoise, + ImagePlaneBackground, + DetectorModePropertiesSetter, + SelectorWheel, ) from scopesim.utils import find_file @@ -197,6 +202,33 @@ def test_ignore_effects_works(self, cmds_with_ignore): class TestTraceDetectorCoordinates: + def test_explicit_detector_mode_resolves_its_binning_before_mapping(self): + train, _ = spectral_geometry_train() + before = train.trace_detector_coordinates() + train.cmds = UserCommands(properties={'!OBS.detector_readout_mode': 'slow'}) + setter = DetectorModePropertiesSetter(mode_properties={'slow': {'!DET.binx': 2}}, cmds=train.cmds) + train.optics_manager.detector_array_effects = [setter] + train.optics_manager.detector_effects = [OnChipBinning(binx='!DET.binx', biny=1, cmds=train.cmds)] + after = train.trace_detector_coordinates() + np.testing.assert_allclose(after['detector_x'], (before['detector_x'] + .5*u.pixel)/2 - .5*u.pixel) + np.testing.assert_array_equal(after['detector_y'], before['detector_y']) + + @pytest.mark.parametrize('selected', [False, True]) + def test_additive_background_preserves_trace_coordinates(self, selected): + train, _ = spectral_geometry_train() + before = train.trace_detector_coordinates() + background = SelectorWheel(selector_key='image_plane_id', wheel=[{'selector_value': [0], 'effect_class': 'ImagePlaneBackground', 'effect_kwargs': {'value': 1.}}]) if selected else ImagePlaneBackground(value=1.) + train.optics_manager.image_plane_effects = [background] + after = train.trace_detector_coordinates() + np.testing.assert_array_equal(before['detector_x'], after['detector_x']) + np.testing.assert_array_equal(before['detector_y'], after['detector_y']) + + def test_unsupported_image_plane_transform_is_not_silently_ignored(self): + train, _ = spectral_geometry_train() + train.optics_manager.image_plane_effects = [Effect()] + with pytest.raises(NotImplementedError, match='additive'): + train.trace_detector_coordinates() + def test_maps_live_trace_to_zero_origin_detector_pixels(self): train, _ = spectral_geometry_train() @@ -258,14 +290,19 @@ def test_rejects_multiple_detectors_on_one_image_plane(self): with pytest.raises(NotImplementedError, match="multiple detectors"): train.trace_detector_coordinates() - def test_rejects_post_extraction_binning(self): + @pytest.mark.parametrize('binning_class', [UnequalBinnedImage, OnChipBinning]) + def test_binned_coordinates_match_before_and_after_readout_geometry(self, binning_class): train, _ = spectral_geometry_train() - train.optics_manager.detector_effects = [ - UnequalBinnedImage(binx=2, biny=1), - ] - - with pytest.raises(NotImplementedError, match="binning"): - train.trace_detector_coordinates() + native = train.trace_detector_coordinates() + binning = binning_class(binx=2, biny=1) + train.optics_manager.detector_effects = [binning] + before = train.trace_detector_coordinates() + np.testing.assert_allclose(before['detector_x'], (native['detector_x'] + .5*u.pixel)/2 - .5*u.pixel) + np.testing.assert_allclose(before['detector_y'], native['detector_y']) + binning.apply_to(train.detector_managers[0][0]) + after = train.trace_detector_coordinates() + np.testing.assert_array_equal(before['detector_x'], after['detector_x']) + np.testing.assert_array_equal(before['detector_y'], after['detector_y']) @pytest.mark.slow @@ -378,6 +415,28 @@ def test_works_with_a_pointer_to_fits_imagehdu(self, cmds): @pytest.mark.usefixtures("patch_mock_path") class TestReadout: + def test_explicit_readout_keys_and_reset(self, unity_cmds): + opt = OpticalTrain(deepcopy(unity_cmds)) + for seconds in (3600., 300., 3600.): + opt.readout(**{'!OBS.dit': seconds, '!OBS.ndit': 1, '!OBS.dit_blue': seconds / 2}) + assert opt.cmds['!OBS.dit'] == seconds + assert opt.cmds['!OBS.dit_blue'] == seconds / 2 + opt.readout(dit=17., ndit=2) + assert opt.cmds['!OBS.dit'] == 17. + assert opt.cmds['!OBS.ndit'] == 2 + + def test_successive_noise_realizations_can_be_replayed(self, unity_cmds, unity_src): + commands = deepcopy(unity_cmds) + commands['!SIM.random.seed'] = 17 + opt = OpticalTrain(commands) + opt.optics_manager.add_effect(BasicReadoutNoise(noise_std=3., ndit=1, cmds=opt.cmds)) + opt.observe(unity_src) + first = opt.readout()[0][1].data + second = opt.readout()[0][1].data + replay = opt.readout(roid=0)[0][1].data + assert not np.array_equal(first, second) + np.testing.assert_array_equal(first, replay) + def test_readout_works_when_source_observed(self, unity_cmds, unity_src): opt = OpticalTrain(unity_cmds)