[10/29] Model configured avalanche noise and warn at detector exposure limits

scopesim/effects/electronic/__init__.py

@@ -52,8 +52,8 @@ def apply_to(self, obj, **kwargs):
 
 from .electrons import LinearityCurve, ADConversion, InterPixelCapacitance
 from .noise import (Bias, PoorMansHxRGReadoutNoise, BasicReadoutNoise,
-                    PixelResponseNonUniformity, ShotNoise, DarkCurrent, ReadoutNoiseFromTable)
-from .exposure import AutoExposure, ExposureIntegration, ExposureOutput
+                    PixelResponseNonUniformity, ShotNoise, AvalancheShotNoise, DarkCurrent, ReadoutNoiseFromTable)
+from .exposure import AutoExposure, ExposureIntegration, ExposureOutput, ExposureLimit
 from .pixels import ReferencePixelBorder, BinnedImage, UnequalBinnedImage, OnChipBinning
 from .dmps import DetectorModePropertiesSetter
 

scopesim/effects/electronic/exposure.py

@@ -2,6 +2,7 @@
 """Exposure actions."""
 
 from typing import ClassVar
+from warnings import warn
 
 import numpy as np
 
@@ -293,3 +294,32 @@ def apply_to(self, obj, **kwargs):
         obj.data = obj.data * dit * ndit
 
         return obj
+
+
+
+class ExposureLimit(Effect):
+    """Warn about expected charge per integration without changing the exposure.
+
+    Apply after integration, dark current and on-chip binning, before noise.
+    The configured usable well refers to a charge packet at this readout stage.
+    Software binning is deliberately later and does not change this limit.
+    """
+
+    required_keys = {"usable_well", "ndit"}
+    z_order = (819,)
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        check_keys(self.meta, self.required_keys, action="error")
+
+    def apply_to(self, obj, **kwargs):
+        if isinstance(obj, Detector):
+            well = from_currsys(self.meta["usable_well"], self.cmds)
+            ndit = from_currsys(self.meta["ndit"], self.cmds)
+            peak = obj.data.max() / ndit
+            obj.header["PEAKDIT"] = (peak, "Expected peak charge per integration [electron]")
+            obj.meta["WELLDEP"] = well
+            obj.header["WELLDEP"] = (well, "Usable charge per integration [input electron]")
+            if peak > well:
+                warn(f"Detector {obj.det_id}: peak {peak:.1f} e per integration exceeds usable well {well:.1f} e; DIT/NDIT are unchanged.", UserWarning, stacklevel=2)
+        return obj

scopesim/effects/electronic/noise.py

@@ -141,10 +141,9 @@ def create_rng(self, det_id: int | None = None) -> np.random.Generator:
             New Generator instance.
 
         """
-        if det_id is not None:
-            seed = [det_id, *self.random_seed]
-        else:
-            seed = self.random_seed
+        seed = self.random_seed
+        if det_id is not None and seed is not None:
+            seed = [det_id, *seed]
         return np.random.default_rng(seed)
 
     def plot(self, det):
@@ -194,7 +193,7 @@ def _apply_to_det(self, det: Detector) -> None:
         logger.debug("Apply %s to %s", self.display_name, det)
 
         ndit = from_currsys(self.meta["ndit"], self.cmds)
-        det.data = self(det.data, ndit)
+        det.data = self(det.data, ndit, det_id=det.det_id)
 
 
 class PoorMansHxRGReadoutNoise(BasicReadoutNoise):
@@ -427,6 +426,9 @@ class ShotNoise(RandomEffect):
     """
     z_order: ClassVar[tuple[int, ...]] = (820,)
 
+    def _apply_to_det(self, det):
+        det.data = self(det.data, det_id=det.det_id)
+
     def __call__(self, data: ArrayLike, det_id: int | None = None) -> NDArray:
         rng = self.create_rng(det_id)
 
@@ -449,3 +451,37 @@ def __call__(self, data: ArrayLike, det_id: int | None = None) -> NDArray:
         )
 
         return data
+
+
+
+class AvalancheShotNoise(ShotNoise):
+    """Input-referred photon multiplication noise, with independent dark charge.
+
+    Conditional on N detected photons, a gamma variate of mean N and variance
+    N(F-1) gives total photon variance F*mu. F is a variance multiplier, not
+    an RMS factor. The adopted input-referred dark term remains Poisson.
+    Apply after ExposureIntegration and DarkCurrent, before read noise.
+    """
+
+    required_keys = {"excess_noise_factor", "dark_current", "dit", "ndit"}
+
+    @property
+    def excess_noise_factor(self):
+        return from_currsys(self.meta["excess_noise_factor"], self.cmds)
+
+    def __call__(self, data, dark_charge=0., det_id=None):
+        factor = self.excess_noise_factor
+        if factor < 1:
+            raise ValueError("The avalanche variance multiplier must be at least one.")
+        if factor == 1:
+            return super().__call__(data, det_id=det_id)
+        rng = self.create_rng(det_id)
+        photons = rng.poisson(data - dark_charge)
+        signal = rng.gamma(photons / (factor - 1), factor - 1)
+        return signal + rng.poisson(dark_charge, size=data.shape)
+
+    def _apply_to_det(self, det):
+        dark = from_currsys(self.meta["dark_current"], self.cmds)
+        dit = from_currsys(self.meta["dit"], self.cmds)
+        ndit = from_currsys(self.meta["ndit"], self.cmds)
+        det.data = self(det.data, dark * dit * ndit, det_id=det.det_id)

scopesim/effects/illumination.py

@@ -212,6 +212,7 @@ def _effect_from_spec(spec, cmds=None):
 
 
 def _detector_qe_from_spec(spec, filename=None, cmds=None):
+    spec = from_currsys(spec, cmds)
     if spec is not None:
         if isinstance(spec, str):
             return SpectralQuantumEfficiency(filename=spec, cmds=cmds)

scopesim/effects/selector_wheel.py

@@ -13,12 +13,13 @@
 import importlib
 from numbers import Integral
 
-from ..utils import (check_keys, get_logger, real_colname)
+from ..utils import from_currsys, check_keys, get_logger, real_colname
 from .effects import Effect
 from ..optics.fov_volume_list import FovVolumeList
 from ..optics.fov import FieldOfView
 from ..optics.image_plane import ImagePlane
 from ..detector.detector import Detector
+from ..source.source import Source
 
 logger = get_logger(__name__)
 
@@ -62,8 +63,10 @@ def __init__(self, **kwargs):
         self.wheel_effects = {}
         for wheel_entry in self.meta["wheel"]:
             selector_value = wheel_entry["selector_value"] # can be single value or list of values
-            effect_class_name = wheel_entry["effect_class"]
+            effect_class_name = from_currsys(wheel_entry["effect_class"], self.cmds)
             effect_kwargs = wheel_entry.get("effect_kwargs", {})
+            if isinstance(effect_kwargs, str):
+                effect_kwargs = from_currsys(effect_kwargs, self.cmds)
 
             # Dynamically import the effect class
             effect_module = importlib.import_module("scopesim.effects")
@@ -78,8 +81,14 @@ def __init__(self, **kwargs):
         self.z_order = self._resolve_z_order()
 
 
-    def apply_to(self, obj, **kwargs):
-        """Based on the selector_key's value in obj.meta, apply the corresponding effect from the wheel."""
+    def apply_to(self, obj, fov=None, **kwargs):
+        """Select from obj.meta, or fov.meta when applying branch effects to a Source."""
+        if isinstance(obj, Source) and fov is not None:
+            selector_value = fov.meta[self.meta["selector_key"]]
+            effect_to_apply = self.get_effect(selector_value)
+            if effect_to_apply is not None:
+                obj = effect_to_apply.apply_to(obj, fov=fov, **kwargs)
+
         if isinstance(obj, FieldOfView):
             if self.meta['selector_key'] not in obj.meta.keys():
                 raise ValueError(f"Selector key {self.meta['selector_key']} not found in FieldOfView meta.")

scopesim/effects/ter_curves.py

@@ -4,6 +4,7 @@
 import warnings
 from typing import ClassVar
 from collections.abc import Collection, Iterable
+from numbers import Integral
 import numpy as np
 from numpy.typing import ArrayLike, NDArray
 from scipy import integrate
@@ -1564,13 +1565,17 @@ class DichroicTree(Effect):
     ```
 
     This effect creates a volume of FoVs, one for each row in the dichroic tree table.
+    For imaging, pass a Source with its fov context to apply the branch throughput before rendering.
     """
     required_keys = {"tree_filename", "dichroic_names", "filename_format"}
     z_order: ClassVar[tuple[int, ...]] = (117, 217, 617)
 
     def __init__(self, cmds=None, **kwargs):
         check_keys(kwargs, self.required_keys, action="error")
         super().__init__(cmds=cmds, **kwargs)
+        if "z_order" in self.meta:
+            z_order = self.meta["z_order"]
+            self.z_order = (int(z_order),) if isinstance(z_order, Integral) else tuple(z_order)
 
         # Load all the dichroic TER curves
         self.dichroics = {}
@@ -1587,7 +1592,7 @@ def __init__(self, cmds=None, **kwargs):
         tree = DataContainer(filename=tree_path)
         self.table = tree.table
 
-    def apply_to(self, obj, **kwargs):
+    def apply_to(self, obj, fov=None, **kwargs):
         """
         To the object, loop through each row of the table and apply the dichroics (defined by columns)
         with the corresponding action defined in the table cell.
@@ -1625,9 +1630,11 @@ def apply_to(self, obj, **kwargs):
             obj.volumes = new_vols
 
         # 2. During observe
-        if isinstance(obj, FieldOfView):
+        if isinstance(obj, FieldOfView) or (isinstance(obj, Source) and fov is not None):
             logger.debug("Executing %s, observe", self.meta['name'])
-            obj_id = obj.meta.get(id_colname, None)
+            branch_fov = obj if isinstance(obj, FieldOfView) else fov
+            obj_id = branch_fov.meta[id_colname]
+            setattr(branch_fov, id_colname, obj_id)
             row = self.table[self.table[id_colname] == obj_id]
 
             if len(row) != 1:

scopesim/optics/optical_train.py

@@ -38,8 +38,14 @@
 N_PROCESSES = mp.cpu_count() - 1
 USE_MULTIPROCESSING = False
 
-def extract_source(fov, source):
-    fov.extract_from(source)
+def extract_source(fov, source, source_effects=(), spatial_header=None, _spectral_wavesets=None):
+    if source_effects:
+        source = source.make_copy()
+        for effect in source_effects:
+            source = effect.apply_to(source, fov=fov)
+        # Branch filtering changes the spectrum and its native wavelength grid.
+        _spectral_wavesets = None
+    fov.extract_from(source, spatial_header=spatial_header, _spectral_wavesets=_spectral_wavesets)
     return fov
 
 def view_fov(fov, hdu_type):
@@ -535,7 +541,15 @@ def observe(self, orig_source=None, update=True, **kwargs):
         source = self.prepare_source(source)
 
         # [1D - transmission curves]
-        for effect in self.optics_manager.source_effects:
+        source_effects = self.optics_manager.source_effects
+        branch_effects = []
+        if not self.fov_manager.is_spectroscope:
+            for ii, effect in enumerate(source_effects):
+                if isinstance(effect, (DichroicTree, SelectorWheel)):
+                    # Attenuation and emission downstream of a split belong to each branch.
+                    source_effects, branch_effects = source_effects[:ii], source_effects[ii:]
+                    break
+        for effect in source_effects:
             source = effect.apply_to(source)
 
         # This source is a private copy, after all source effects. FOV extraction
@@ -558,7 +572,7 @@ def observe(self, orig_source=None, update=True, **kwargs):
 
             with mp.Pool(processes=N_PROCESSES) as pool:
                 fovs = pool.starmap(extract_source,
-                                    zip(fovs, [source] * len(fovs)))
+                                    zip(fovs, [source] * len(fovs), [branch_effects] * len(fovs)))
 
             with mp.Pool(processes=N_PROCESSES) as pool:
                 fovs = pool.starmap(view_fov,
@@ -614,9 +628,9 @@ def observe(self, orig_source=None, update=True, **kwargs):
                     spatial_header = foveffs[psf_index].padded_header(fov, point_sources=point_sources)
                     for effect in spatial_shifts:
                         spatial_header = effect._padded_header(spatial_header, effect.get_shifts(fov))
-                    fov.extract_from(source, spatial_header=spatial_header, _spectral_wavesets=spectral_wavesets)
+                    extract_source(fov, source, branch_effects, spatial_header=spatial_header, _spectral_wavesets=spectral_wavesets)
                 else:
-                    fov.extract_from(source, _spectral_wavesets=spectral_wavesets)
+                    extract_source(fov, source, branch_effects, _spectral_wavesets=spectral_wavesets)
                     fov.view()
                 nobar = len(foveffs) <= 1
                 for effect in tqdm(foveffs, disable=nobar,

scopesim/tests/tests_effects/test_MoffatPSF.py

@@ -1,4 +1,5 @@
 import logging
+from unittest.mock import Mock
 
 import numpy as np
 import pytest
@@ -251,3 +252,33 @@ def test_ao_cache_tracks_resolved_path_contents_and_alpha(tmp_path):
     np.testing.assert_allclose(psf.fwhm(wave).value, .2)
     psf.meta['alpha'] = 3.
     assert psf.alpha == 3.
+
+
+@pytest.mark.parametrize("psf_class", [MoffatPSF, AOEnhanceablePSF])
+def test_observation_configuration_is_only_required_when_evaluating(psf_class, monkeypatch):
+    observation_info = Mock(side_effect=KeyError("!OBS.airmass"))
+    monkeypatch.setattr("scopesim.effects.psfs.analytical.get_observation_info_from_cmds", observation_info)
+    kwargs = {"alpha": 4.765, "fwhm": {"seeing": 0.7, "seeing_unit": "arcsec", "pivot_wave": 500, "pivot_wave_unit": "nm"}}
+    if psf_class is AOEnhanceablePSF:
+        kwargs.update(ao_table="unused.dat", enable_ao=False, is_absolute=False)
+    psf = psf_class(**kwargs)
+    observation_info.assert_not_called()
+    with pytest.raises(KeyError, match="airmass"):
+        psf.fwhm(np.array([0.5]) * u.um)
+
+
+
+def test_ao_toggle_preserves_natural_and_absolute_widths(tmp_path):
+    filename = tmp_path / "ao.ecsv"
+    table = Table({"wavelength": np.array([0.5, 1.0, 1.5, 2.0]) * u.um, "fwhm": [0.5, 0.5, 0.5, 0.5]})
+    table.write(filename, format="ascii.ecsv")
+    psf = AOEnhanceablePSF(alpha=4.765, fwhm=0.7, ao_table=str(filename), enable_ao=False, is_absolute=False)
+    wavelengths = np.array([1.0]) * u.um
+    np.testing.assert_allclose(psf.fwhm(wavelengths).to_value(u.arcsec), [0.7])
+    psf.meta["enable_ao"] = True
+    np.testing.assert_allclose(psf.fwhm(wavelengths).to_value(u.arcsec), [0.35])
+    psf.meta["is_absolute"] = True
+    np.testing.assert_allclose(psf.fwhm(wavelengths).to_value(u.arcsec), [0.5])
+    psf.meta["fwhm"] = 0.9
+    psf.meta["enable_ao"] = False
+    np.testing.assert_allclose(psf.fwhm(wavelengths).to_value(u.arcsec), [0.9])