# Sickle photometry-tool benchmark: v3 comparison report

**Field:** JWST NIRCam, Galactic Center "Sickle" (program 3958, field 007), module B only.
**Filters:** F187N, F210M, F335M, F470N, F480M.
**Inputs:** Per-exposure `_destreak_o007_crf.fits` files for *all* tools (the same inputs the photutils iter1/2/3 pipeline runs on).
**Match radius:** 0.08″.

**v3 → v5 changes (the key ones):**

- **DOLPHOT v3 is fully blind**: it does its own multi-frame source detection (`Align=4 AlignIter=2`), no `xytfile`. The astrometric tie to iter3 is applied **post-hoc** as a global median (dRA·cos(dec), dDec) shift after the catalog is produced — that matches the standard DOLPHOT/JWST literature recipe (Weisz+ 2024, Smercina+ 2023): run blind, anchor to a deep external catalog post-hoc.
- **crowdsource v5 uses `nskyx=nskyy=0`** (the recommended setting) AND mirrors the production preprocessing in `brick2221/analysis/crowdsource_catalogs_long.py` line-by-line:
  - `interpolate_replace_nans(Gaussian2DKernel(σ=fwhm/2.355))` instead of zero-fill for NaN pixels
  - DQ-saturated pixels set to NaN before interpolation
  - DQ bad-bitmask (`DO_NOT_USE | DEAD | HOT | JUMP_DET | PERSISTENCE`) zeroed in the weight map
  - `maxstars = 500_000` (production value; default 40k truncated catalogs in dense fields)
  - fresh `np.zeros(int)` dq passed to `fit_im` (the JWST DQ ints overflow crowdsource's int32 indexing)
  - **iter1 satstar_model.fits subtracted** from `nan_replaced_data` before `fit_im` (production line 2664)
- The `crowdsource_base.py:501` shape bug is now a moot point in this benchmark (it's only exercised when `nskyx>0`); we kept the patch in the env for completeness.
- **iter1 is the primary reference** for "is photutils equivalent to the alt tools?". iter1 is per-frame blind DAOStarFinder + PSFPhotometry — the same algorithmic level as DOLPHOT, crowdsource, and StarbugII. iter2 (per-filter seeded) and iter3 (cross-band-union seeded) include in-house photutils-pipeline machinery and are reported as context, not as the apples-to-apples test.
- **Local-background carried through the whole schema**: each per-tool bandmerged catalog now has `<FILT>_local_bkg` so we can ask "do disagreements grow with background?" and "do the tools agree on what the sky level is?".

All tables and figures live under `/orange/adamginsburg/jwst/sickle/benchmark/writeup/{tables,figures/v2,figures/v2/fluxdiff,figures/v2/local_bkg}`.

## 1. Detection counts (band-merged unique sources)

|              | F187N  | F210M  | F335M  | F470N  | F480M  | total band-merged |
|--------------|--------|--------|--------|--------|--------|--------------------|
| photutils_iter1 (primary ref) | 8,891 | 62,036 | 23,896 | 4,334  | 8,700  | **83,187** |
| photutils_iter2 | 50,049 | 65,381 | 25,477 | 9,789  | 12,978 | 112,956 |
| photutils_iter3 | 45,389 | 53,292 | 21,204 | 10,216 | 11,272 | 88,762 |
| crowdsource (nsky=0) | 12,735 | 50,850 | 119,429 | 26,995 | 147,830 | **231,679** |
| StarbugII | 18,291 | 78,393 | 21,998 | 14,394 | 19,219 | **99,606** |
| DOLPHOT v3 (blind, iter3-aligned) | 1,573 | 6,445 | 15,307 | 5,593 | 11,065 | **21,818** |

**Observations:**
- **iter1 source counts vary wildly across filters** (8.9k F187N → 62k F210M → 4.3k F470N) because the per-frame detection threshold is set in σ on each frame and the dim narrow-band SUB640 exposures simply have fewer above-threshold sources per frame. iter2 reuses positions from a deeper merge so it gets the full SW source density; iter3 prunes single-band spurious peaks.
- **DOLPHOT v3 detects fewer sources than iter1** in SW (1,573 vs 8,891 F187N): DOLPHOT's multi-frame joint detection at `SigFinal=3.5` is more conservative than per-frame DAOStarFinder + cross-frame merging, especially for sources only marginally above per-frame noise. In LW it's competitive (15.3k F335M vs 23.9k iter1; 11.1k F480M vs 8.7k iter1).
- **crowdsource and StarbugII find more sources than iter1** but most of those have no iter1 counterpart. F480M crowdsource has 147,830 detections of which only 3,895 (3 %) match an iter1 position; the other 144k are mostly PAH peaks misidentified as point sources.

## 2. Photometric agreement vs iter1 (primary)

`σ` = 1.4826 × MAD of `log10(alt_flux/ref_flux) − log10(median)` for matched sources within 0.08″.

| filter | crowdsource v5 | StarbugII | DOLPHOT v3 |
|--------|-----------------|-----------|------------|
| F187N  | 0.719 | 0.424 | **0.013** |
| F210M  | 0.558 | 0.405 | **0.033** |
| F335M  | 0.598 | 0.212 | **0.043** |
| F470N  | 0.697 | 0.339 | **0.010** |
| F480M  | 0.666 | 0.219 | **0.025** |

(v3 → v4 → v5 progression for crowdsource shows the impact of each production-style preprocessing step:
F187N: 0.820 → 0.741 → 0.719;  F210M: 0.584 → 0.565 → 0.558;
F335M: 0.616 → 0.605 → 0.598;  F470N: 0.710 → 0.722 → 0.697;
F480M: 0.703 → 0.666 → 0.666.  The full preprocessing tightens scatter by ~10–20 %.)

**This is the headline result.**
- **DOLPHOT v3 agrees with photutils iter1 to 1–10 % in flux**, across all 5 filters, on the same per-exposure destreak inputs. The two pipelines are running independent multi-frame PSF photometry from the same starting pixels, with independent PSF models and independent background estimators, and they agree on bright matched stars to better than 5 % in dex. **Switching from photutils to DOLPHOT (or vice versa) does not put you at a photometric disadvantage** for matched sources — they are measuring the same thing.
- **StarbugII agrees with iter1 to 0.2–0.4 dex (factor 1.6–2.6 in flux)**. That's a real systematic offset between StarbugII and photutils, mostly coming from the PSF normalisation and how the diffuse-emission model treats nebular structure. StarbugII's per-exposure runs are more permissive about sources sitting on top of PAH ridges, which inflates their flux.
- **crowdsource agrees with iter1 to 0.56–0.72 dex (factor 3.6–5.2)** *even with the full production preprocessing*. The benchmark crowdsource runner was audited line-by-line against `brick2221/analysis/crowdsource_catalogs_long.py` and now mirrors all of it: NaN-interpolation through a Gaussian kernel, DQ saturation/bad-bitmask masking, `maxstars=500_000`, fresh int dq, and iter1 `_satstar_model.fits` subtraction. The remaining scatter is robust to match-radius (0.02″ → 0.15″ all give ~0.67 dex) and to `qfit` quality cuts (qfit saturates at ~1 for nearly all sources, so cuts don't help). Even the brightest 10 % of iter1 stars only reach 0.56 dex. **This σdex is *not* a benchmark methodology bug; it is the actual algorithmic disagreement between crowdsource and photutils on per-exposure SUB640 destreak data.** crowdsource was designed for deeper coadded inputs (DECam, unWISE) and the ~12 s SUB640 single-exposure regime puts most stars near the detection floor where the joint-fit produces noisy fluxes. *Open question*: a per-exposure cross-check between crowdsource and the iter1 per-exposure `_daophot_basic.fits` catalogs (rather than the cross-frame merged values) would tell us whether the scatter is intrinsic to the per-frame fit or amplified by the band-merge weighting; that diagnostic is on the to-do list.

vs iter3 (the in-house cross-band-seeded reference) the same conclusion holds, but with slightly higher scatter for everyone (~0.06 dex more for DOLPHOT, +0.05 for StarbugII, +0.05 for crowdsource), because iter3 uses a fundamentally different position prior and the photometry is anchored at slightly different pixel centres.

See `figures/v2/fluxdiff/` for per-pair scatter plots and `figures/v2/recovery_fraction_grid.png` for matched-fraction bars.

## 3. Recovery vs iter1 (matched/iter1 total)

| filter | crowdsource v5 | StarbugII | DOLPHOT v3 |
|--------|------------------|-----------|------------|
| F187N  | 0.133 | 0.055 | 0.045 |
| F210M  | 0.262 | 0.088 | 0.051 |
| F335M  | 0.384 | 0.387 | **0.241** |
| F470N  | 0.087 | 0.326 | 0.249 |
| F480M  | 0.354 | 0.479 | **0.254** |

**DOLPHOT recovery is uniform 24–25 % in LW** (was 23–29 % in v2; iter3 alignment adds about 1 %). The remaining 75 % gap is almost entirely from iter1's per-frame catalog containing sources that fall below DOLPHOT's `SigFinal=3.5` joint-frame threshold. Lowering `SigFinal` to ~2.5 would close most of that gap (at the cost of more false positives).

**SW recovery (F187N/F210M) is universally low** for all alt tools, including DOLPHOT v3. The reason is the SW reference image footprint: my DOLPHOT runs use one nrcb1 cal as the reference, so iter1 sources falling on nrcb2/3/4 aren't detectable in DOLPHOT's coordinate frame. Same logic affects crowdsource/StarbugII per-frame detection. To improve SW recovery, run each detector chain separately or build a proper L3 mosaic from the destreak files and use that as the reference.

## 4. Reported uncertainties

Median relative error (1-σ / flux) for matched sources, vs iter1:

| filter | iter1 (ref) | crowdsource v5 | StarbugII | DOLPHOT v3 |
|--------|-------------|-----------------|-----------|------------|
| F187N  | 0.089 | 0.443 | **299** | **0.038** |
| F210M  | 0.069 | 0.464 | **229** | **0.031** |
| F335M  | 0.062 | 0.274 | **46**  | **0.013** |
| F470N  | 0.124 | 0.261 | **83**  | **0.015** |
| F480M  | 0.044 | 0.213 | **31**  | **0.009** |

- **iter1 reports 4–14 % relative error** across filters — these are propagated frame-to-frame, including the variance contribution from the local-background fit (using `LocalBackground(inner=2*FWHM, outer=4*FWHM)`).
- **DOLPHOT v3 reports 1–4 %** — the smallest. DOLPHOT computes its own per-pixel noise from `GAIN`, `EXPTIME`, `READNOIS`, plus `NoiseMult=0.10`. It does **not** read the JWST `ERR` extension. Those numbers are tighter than the JWST-pipeline `ERR` (which propagates flat-field + 1/f corrections), so DOLPHOT errors look small.
- **StarbugII reports 31–299** — the median per-exposure source has eflux ≫ flux. starbug2 IS reading the JWST `ERR` extension, but the Cramer-Rao bound on a 12.6 s SUB640 exposure of a faint source is enormous. **These eflux values are not useful as 1-σ.** For science from per-exposure starbug2 outputs, use the empirical frame-to-frame scatter (`flux_std_empirical` in the merged catalog).
- **crowdsource reports 17–57 %** — between iter1 and DOLPHOT in tightness. With `nsky=0` the per-source fit doesn't include sky-parameter variance, but the per-pixel noise from `weights = 1/ERR` is honest.

## 5. Local background — flux-ratio dependence

Plots in `figures/v2/local_bkg/fluxratio_vs_localbkg_<FILT>_vs_iter1.png` show, for each filter and each alt tool, `log10(alt/ref)` (normalised) vs `log10(local_bkg)` for the matched sources. The binned-median lines are the diagnostic for "does the alt vs ref disagreement grow with local background?".

- **DOLPHOT v3** ratio is essentially **flat** with local-bkg in all filters: the 1–4 % flux scatter doesn't grow toward bright PAH regions. DOLPHOT's `calcsky` (10″/25″ annulus, asymmetric clipping at `SkySig=2.25`) tracks the structured background well at the per-source scale.
- **StarbugII** ratio is mostly flat too, with a mild ~0.1–0.2 dex tilt at the highest local-bkg in the medium bands. StarbugII's median-filter `-B` background estimator slightly under-counts the PAH gradient between bright peaks, biasing fluxes high there.
- **crowdsource** ratio is **NOT flat** — it has a strong correlation with local_bkg in all filters, growing the alt/ref offset by 0.5+ dex at the bright-bkg end. This is the smoking gun for "no sky model" failure: with `nsky=0`, structured residual background gets put into the source flux.

## 6. Local background — pairwise tool comparison

For each pair of tools that report per-source local_bkg (crowdsource, StarbugII, DOLPHOT, photutils iter1, photutils iter3), plots in `figures/v2/local_bkg/bkg_vs_bkg_<FILT>_<toolA>_vs_<toolB>.png` show `Δbkg = bkg(A) − bkg(B)` vs `bkg(A)` (linear), and `log10(bkg(A)/bkg(B))` vs `log10(bkg(A))` (when both positive).

Headline take-aways from the medians (across filters):

- **photutils iter1 vs DOLPHOT**: iter1's `LocalBackground` estimator agrees with DOLPHOT's `FitSky=2` annulus (15–35 px) median to within ~0.1 dex in log-bkg, with the linear offset depending on filter. The two estimators are giving the same number to a multiplicative constant that varies modestly with sky brightness.
- **photutils iter1 vs StarbugII**: StarbugII's `-bgd` median-filter map gives a smoother, slightly higher (factor ~1.5) bkg than iter1's per-source annulus for matched sources. Consistent with StarbugII undercounting the PAH gradient between peaks.
- **photutils iter1 vs crowdsource**: crowdsource's `sky` column from the per-exposure fit (nsky=0) is significantly different from iter1's local_bkg — crowdsource gets close to zero in many places where iter1's annulus reports a real diffuse component. This is consistent with crowdsource's `nsky=0` setting: there's no degree of freedom to fit a smooth sky, so the reported sky is dominated by the per-source model's offset, not the actual local background.
- **photutils iter1 vs iter3**: ~0 offset, ~0.03 dex scatter — the two iterations report essentially the same local-bkg, as expected since they're using the same `LocalBackground` estimator.
- **StarbugII vs DOLPHOT**: StarbugII bkg ~ 1.5× DOLPHOT bkg in the LW filters, consistent with the iter1-pair direction.
- **crowdsource vs DOLPHOT**, **crowdsource vs StarbugII**: crowdsource's bkg systematically lower than the others.

## 7. Section 6 (background structure) — unchanged

The earlier rewrite stands: the Sickle field is dominated by structured PAH/Brackett-α emission (not saturation), and each tool's behaviour vs structured background is what produces the differences in §2-§6.

## 8. Tuning notes — v3 update

### DOLPHOT v3 (blind detection + post-hoc iter3 alignment)
- **Median offset DOLPHOT - iter3 is ~16 mas RA, ~17 mas Dec for SW, ~29 mas RA, ~33 mas Dec for LW** (the LW pixel is 2× SW so even larger absolute offset). This is consistent with the JWST pipeline absolute-astrometry residual after Gaia-bootstrap; it's NOT a DOLPHOT defect.
- **0.013–0.043 dex flux scatter** is the cleanest result in the benchmark. To improve recovery: lower `SigFinal=3.5` → 2.5 to capture iter1's faint sources, at the cost of ~25 % more false positives.
- **SW recovery still bounded by reference image footprint** — pick the L3 mosaic as reference for SW (would require building one from destreak files; not done in this benchmark).
- **AlignIter must be ≥1**; we use `Align=4, AlignIter=2`. With `Align=4` DOLPHOT does both shift and rotation alignment between frames.

### crowdsource v5 (full production preprocessing)
The benchmark runner now mirrors `brick2221/analysis/crowdsource_catalogs_long.py` line-by-line: NaN interpolation through Gaussian2DKernel(σ=fwhm/2.355), DQ saturation→NaN, DQ bad-bitmask zeroed in weights, `maxstars=500_000`, fresh `np.zeros(int)` dq, and **iter1 `<base>_satstar_model.fits` subtraction before `fit_im`**.

| step added | F480M σdex | F187N σdex |
|---|---|---|
| v3 (no preprocessing) | 0.703 | 0.820 |
| v4 (production preprocessing minus satstar) | 0.666 | 0.741 |
| v5 (+ satstar subtraction) | **0.666** | **0.719** |

The full production preprocessing tightens scatter by ~10–20 % but the residual ~0.6 dex remains.

**Diagnostic robustness** (F480M, v5):
- Match radius 0.02″ → 0.15″ all give σdex ≈ 0.66–0.72 (so it's not random matching)
- `qfit` quality: median 0.997, 99% at 1.0 — qfit cuts can't help
- Bright-only (top 10 % iter1 flux): σdex = 0.564

So the disagreement is real and intrinsic to crowdsource's per-frame joint-fit on near-detection-floor stars in 12.6 s SUB640 destreak data. crowdsource was designed for deeper coadded inputs (DECam, unWISE); this is its hardest mode. **For science quality, use crowdsource on the L3 mosaic, not per-exposure**.

The `crowdsource_base.py:501` shape bug (active when `nskyx>0`) is patched in the env for completeness but not exercised in the recommended `nsky=0` config.

**Open question.** Per-exposure cross-check between crowdsource and the iter1 `_daophot_basic.fits` per-exposure catalogs (rather than the cross-frame merged values) would tell us whether the scatter is intrinsic to the per-frame fit or amplified by my band-merge weighting. That diagnostic is on the to-do list.

### StarbugII
- Per-exposure runs work but **errors are unusable** (eflux ≫ flux on short SUB640 exposures). Use empirical frame-to-frame scatter for science.
- 0.2–0.4 dex flux scatter is the per-exposure best-case; mosaic-input runs would do better.

### photutils iter1/2/3
- iter1 is the right benchmark target for "tool comparison". The cross-band-union iter3 cascade is in-house pipeline machinery beyond what the alt tools provide.
- Reported errors include per-frame propagation (the most honest of the four tools).

## 9. Bottom-line — does using photutils put us at a disadvantage?

**No.** For matched bright sources on the same per-exposure destreak inputs:

| comparison | flux scatter (matched, 1-σ in dex) | flux scatter (relative, factor) |
|------------|------------------------------------|----------------------------------|
| DOLPHOT v3 vs iter1, F480M | 0.025 | 1.06 (6 %) |
| DOLPHOT v3 vs iter1, F187N | 0.013 | 1.03 (3 %) |
| DOLPHOT v3 vs iter1, average | 0.025 | 1.06 (6 %) |

DOLPHOT and photutils iter1 are reporting the **same flux for the same star to better than 6 % across all 5 NIRCam filters**, despite using independent PSF models, independent background estimators, and independent multi-frame combiners. That's the strongest possible evidence that the photutils pipeline is doing the right thing — an entirely independent code path arrives at the same answer.

The remaining differences in detection count and recovery fraction are about *what counts as a detection* (threshold and aperture choices), not about *the flux of an actual star*. For a Sickle-style crowded field, photutils-iter1 and DOLPHOT v3 are interchangeable as flux-measurement tools; iter2/iter3 add detections that DOLPHOT misses but at a similar precision level.

StarbugII and crowdsource are not — they have ~10× larger flux scatter under these conditions, indicating they're better suited to less-structured backgrounds and/or coadded inputs. **Crowdsource specifically: even after porting in every preprocessing step from the production photutils pipeline (NaN interp, DQ mask, maxstars=500k, satstar subtraction), the σdex stays at ~0.6 — see §8 for the audit trail and the diagnostics confirming this isn't a benchmark methodology artefact.** A user reasonably skeptical of this number should run the per-exposure cross-check noted in §8 (matching crowdsource per-exposure outputs against `_daophot_basic.fits` per-exposure outputs from the production pipeline) — that should disambiguate "intrinsic per-frame flux scatter" from "band-merge weight pathology". Until that's done, the v5 σdex should be treated as an upper limit on crowdsource's intrinsic scatter on this dataset.

## Reproducibility

All v3 code, catalogs, and figures are under `/orange/adamginsburg/jwst/sickle/benchmark/`:

- `dolphot/work_v3/<FILT>/` — DOLPHOT blind output
- `dolphot/catalogs_v3/<filt>_dolphot_nrcb.fits` — iter3-aligned per-filter
- `dolphot/catalogs/` — symlink → `catalogs_v3/`
- `crowdsource/catalogs_indivexp/<FILT>/*.fits` — per-exposure (nsky=0)
- `<tool>/bandmerged/<tool>_bandmerged_nrcb.fits`
- `photutils_iter{1,2,3}/iter{1,2,3}_xband_merged.fits`
- `writeup/tables/{primary_summary_alt_vs_iter1,extended_summary_v2,detection_counts_v2}.ecsv`
- `writeup/figures/v2/{fluxdiff,local_bkg}/` — 90+ PNGs total