Fix fractional scatter height with duplicated traces; re-enable scatter for HARPSPOL

CHANGELOG.md

@@ -4,7 +4,8 @@
 
 ### Fixed
 - **Background scatter is now measured on the frame it corrects.** The `scatter` step fits its model on the flat lamp frames (every instrument's `id_scatter` points at a flat), and `ScienceExtraction` subtracted those coefficients from science frames unchanged. Since a master flat is a *sum* of lamp exposures, the model is in the flux units of that stack and generally has no relation to a single science exposure: measured on XSHOOTER nir (10x400 s flats vs a 5 s science frame) the subtraction was 67x the signal at the trace centres, driving 99.97% of in-order pixels negative; on HARPS red (5x40 s vs 600 s) it *under*-subtracted at 0.14x. `NormalizeFlatField` was unaffected, as it corrects the same flat the model was measured on. Consumers now re-estimate the background on the calibrated image they are about to extract, via the new `ScatterModel.refit`, so the model matches that frame's own inter-order level by construction (1.000 on the XSHOOTER frame above, against 67x before). No exposure-time bookkeeping is involved and no settings changed
-- HARPSPOL: the `scatter` step is no longer declared, having been measured to remove most of the stellar flux at every `extraction_height` tried (#38, #39)
+- **A fractional scatter `extraction_height` collapsed to a zero-height aperture wherever traces are stored twice.** `Trace.run` saves `grouped + raw_traces`, so on an instrument with a `fibers` block every trace is present twice, and where a group holds a single fiber the two copies carry identical polynomials. A fractional height is resolved against the median separation between adjacent traces, which is then 0: nothing was masked and the polynomial was fitted to the order flux itself. On HARPSPOL blue the model came out 8.9x the frame's own inter-order floor, identically at every height from 0.2 to 1.0 — that insensitivity being the visible symptom — and extraction gave a median flux of 621.8 with 27.3% of points negative, against 2027.8 uncorrected. Spacing is now measured from distinct positions only, and a list with no distinct positions warns rather than silently using the 10 px fallback. HARPS, UVES, XSHOOTER and LICK_APF have no `fibers` block, so they never duplicate traces and were unaffected
+- HARPSPOL: the `scatter` step is declared again, at `extraction_height: 0.7`. It was removed in #38 on measurements taken with both the pre-refit method and the spacing defect above in play; with each addressed, model/floor is 1.18x on blue and 1.02x on red over a 0.2–1.0 sweep on gamma Equ (HD 201601, 2012-07-16), against 1.33x for HARPS (#38, #39, #40)
 - Scatter `extraction_height` was narrower than the order footprint on two instruments, leaving order wings in the background fit and biasing it high: HARPS shipped 20 px against 77 px order spacing (now `0.9` of the spacing, model/truth 6.1x -> 1.3x) and LICK_APF 18 px with degree 4 (now `0.6` and `scatter_degree: 2`, 0.48x -> 1.2x, the degree drop because 22 px order spacing leaves too few inter-order pixels for a quartic). UVES moved 0.9 -> 1.0 (1.23x -> 1.04x). Measured against the median of pixels further than 0.4x the order spacing from any trace. See `scatter.md`
 
 ### Changed

pyreduce/estimate_background_scatter.py

@@ -120,13 +120,26 @@ def estimate_background_scatter(
     # Compute extraction height in pixels if fractional
     xwd = extraction_height
     if xwd is not None and xwd < 3:
-        # Fraction of order spacing - estimate from trace separation
+        # Fraction of order spacing - estimate from trace separation.
+        # Trace.run stores `grouped + raw_traces`, so on an instrument with a
+        # fibers block every trace is present twice; where a group holds a single
+        # fiber the two copies have identical polynomials. Coincident traces
+        # contribute zero separations, which would drag the median to 0 and
+        # silently collapse the aperture to zero height, masking nothing.
         x_mid = ncol // 2
         y_mids = np.array([np.polyval(t.pos, x_mid) for t in traces])
-        if len(y_mids) > 1:
-            spacing = np.median(np.abs(np.diff(np.sort(y_mids))))
+        separations = np.diff(np.sort(y_mids))
+        separations = separations[separations > 0]
+        if len(separations) > 0:
+            spacing = np.median(separations)
             xwd = xwd * spacing
         else:
+            logger.warning(
+                "Cannot measure order spacing from %d trace(s) with no distinct "
+                "positions; falling back to a %d px aperture for the scatter mask.",
+                len(traces),
+                10,
+            )
             xwd = 10  # fallback
 
     # Method 1: Select all pixels, but those known to be in traces

pyreduce/instruments/HARPSPOL/__init__.py

@@ -59,6 +59,11 @@ def get_expected_values(
                 "night": night,
                 "type": r"(LAMP,LAMP),TUN",
             },
+            "scatter": {
+                "instrument": "HARPS",
+                "night": night,
+                "type": r"(LAMP,LAMP),TUN",
+            },
             "curvature": {
                 "instrument": "HARPS",
                 "night": night,

pyreduce/instruments/HARPSPOL/settings.json

@@ -10,6 +10,9 @@
         "manual": false,
         "border_width": [10, 10, 0, 0]
     },
+    "scatter": {
+        "extraction_height": 0.7
+    },
     "norm_flat": {
         "smooth_slitfunction": 2,
         "extraction_height": 20,

test/test_scatter.py

@@ -148,3 +148,45 @@ def test_masks_traces_whose_valid_columns_are_not_contiguous():
         img, traces, extraction_height=20, scatter_degree=0, border_width=0
     )
     assert np.allclose(coeff[0, 0], 10.0)
+
+
+@pytest.mark.unit
+def test_fractional_height_with_traces_stored_twice():
+    """Coincident traces must not collapse the fractional-height aperture.
+
+    Trace.run stores ``grouped + raw_traces``, so on an instrument with a fibers
+    block every trace is present twice, and where a group holds a single fiber the
+    two copies carry identical polynomials. The zero separations between those
+    copies used to drag the median order spacing to 0, so a fractional
+    extraction_height became a zero-height aperture, nothing was masked, and the
+    polynomial was fitted to the order flux itself.
+    """
+    nrow, ncol = 100, 200
+    background, order_flux = 10.0, 40.0
+    img = np.full((nrow, ncol), background)
+    rows = [25, 50, 75]  # spacing 25, so a 0.8 fraction is a 20 px aperture
+    for r in rows:
+        img[r - 10 : r + 10] = order_flux
+
+    # Order flux this close to the background survives the sigma clip, so an
+    # unmasked fit lands on the area-weighted mean (28.0) rather than the floor.
+    unique = [
+        Trace(m=i, group="upper", pos=np.array([0.0, float(r)]), column_range=(0, ncol))
+        for i, r in enumerate(rows)
+    ]
+    # the same three traces again, as the raw per-fiber copies
+    duplicated = unique + [
+        Trace(m=t.m, group=None, pos=t.pos.copy(), column_range=t.column_range)
+        for t in unique
+    ]
+
+    coeff = estimate_background_scatter(
+        img, duplicated, extraction_height=0.8, scatter_degree=0, border_width=0
+    )
+    assert np.allclose(coeff[0, 0], background)
+
+    # and it agrees with the same fit on the de-duplicated list
+    expected = estimate_background_scatter(
+        img, unique, extraction_height=0.8, scatter_degree=0, border_width=0
+    )
+    assert np.allclose(coeff[0, 0], expected[0, 0])