# Sickle photometry-tool benchmark: comparison report

**Field:** JWST NIRCam, Galactic Center "Sickle" (program 3958, field 007), module B only.
**Filters:** F187N, F210M, F335M, F470N, F480M (5 NRCB).
**Reference:** photutils-based three-iteration cascade ("iter3"), the in-house pipeline at
`brick2221/analysis/crowdsource_catalogs_long.py`, cross-band-union-seeded.
**Other tools tested:** DOLPHOT 2.1 + NIRCam module (Dolphin / Weisz+24), Schlafly's `crowdsourcephoto` 0.5.9, StarbugII 0.7.7 (Nally).
**StarFinder:** skipped on user instruction (IDL, not maintained for JWST).
**Match radius:** 0.08″ (~half the worst-filter FWHM). Pairwise matches against the photutils-iter3 reference; no external ground truth.

All catalogs, residuals, and figures live under `/orange/adamginsburg/jwst/sickle/benchmark/`.

## 1. Data inputs

DOLPHOT consumed the per-exposure `_cal.fits` files and ran a single multi-frame fit per filter (24 LW exposures or 96 SW exposures × 4 detectors). Crowdsource and StarbugII both ran on the level-3 NRCB drizzled mosaic (`{FILT}/pipeline/jw03958-o007_t001_nircam_clear-{filt}-nrcb_i2d.fits`). The photutils iter3 reference is built from the same per-exposure `_destreak_o007_crf.fits` images that the iter1/iter2 pipeline already used, then cross-filter merged via FoF clustering with a 0.041″ link.

Per-filter F187N and F210M cover four short-wave detectors; F335M, F470N, F480M cover only NRCBLONG (sub-array SUB640 in this dataset, 640×640 pixels).

## 2. Detection counts (band-merged, unique sources)

|              | F187N  | F210M  | F335M  | F470N  | F480M  | total band-merged |
|--------------|--------|--------|--------|--------|--------|--------------------|
| photutils_iter3 (ref) | 45,389 | 53,292 | 21,204 | 10,216 | 11,272 | **88,762** |
| crowdsource  | 27,145 | 59,434 | 37,645 | 13,059 | 36,983 | 119,136 |
| StarbugII    | 20,069 | 60,856 | 12,362 |  9,551 | 10,830 | 65,338 |
| DOLPHOT      |  1,692 |  6,386 | 14,826 |  5,792 | 10,613 | 21,460 |

**Take-aways:**
- **Crowdsource** finds the most sources overall, but most of them have no photutils counterpart (≈25-35k tool-only per filter — it is much more permissive in the diffuse PAH-dominated background).
- **DOLPHOT** is by far the most conservative; it deletes ~80 % of F187N sources via its sharpness/chi cuts and PSF S/N threshold (`SigFinal=3.5`). This is a **feature**: in a Galactic Center field most "sources" below 5σ are background structure, not stars, and DOLPHOT's filtering keeps them out. The flip side is that 88 % of the photutils F187N catalog has no DOLPHOT counterpart — DOLPHOT is missing real faint stars (or photutils is keeping spurious ones).
- **StarbugII** and **photutils** sit in the middle. StarbugII actually exceeds photutils in F210M (60,856 vs 53,292) — the strongest agreement of the three tools.

## 3. Recovery fractions (matched / photutils total)

See `figures/recovery_fraction.png`. Numerical values (`tables/extended_summary.ecsv`):

| filter | crowdsource | StarbugII | DOLPHOT |
|--------|-------------|-----------|---------|
| F187N  | 0.037 | 0.34 | 0.024 |
| F210M  | 0.074 | 0.75 | 0.048 |
| F335M  | 0.042 | 0.46 | 0.245 |
| F470N  | 0.017 | 0.54 | 0.20  |
| F480M  | 0.044 | 0.58 | 0.25  |

**Take-aways:**
- **StarbugII tracks photutils best** — recovers 34-75 % of the iter3 sources, depending on filter, with consistent matching even in the medium-band filters where blending is severe.
- **DOLPHOT's recovery is terrible in SW (F187N/F210M)** — only 2-5 % matched. This points at WCS / astrometric registration disagreement at SW-pixel scale, since DOLPHOT itself detected 6,386 sources in F210M but only 2,550 of them landed within 0.08″ of an iter3 source. The DOLPHOT internal alignment uses the first cal as the chip reference and aligns every other frame to that; the photutils pipeline aligns to a Gaia/VVV refcat. Probable next investigation: re-run DOLPHOT with `UseWCS=2 + Align=4` (currently set) but with a known-good `xytfile` from the iter3 catalog.
- **Crowdsource recovery is the worst** — 1.7-7.4 %. Most of crowdsource's sources are **not** at iter3 positions; they cluster on PAH structure that the SimplePSF-wrapped grid PSF misfits. This is the main story for crowdsource in this field (see §4 below).

## 4. Photometric agreement (matched sources only)

Each tool reports flux in its native unit (DOLPHOT: count rate; StarbugII / photutils: MJy/sr × pixel area; crowdsource: same as photutils since both use stpsf grid PSFs). The relevant quantity is the *spread* in the log-flux ratio after re-normalising by the median offset, expressed in dex:

| tool        | F187N σ (dex) | F210M σ | F335M σ | F470N σ | F480M σ |
|-------------|---------------|---------|---------|---------|---------|
| crowdsource | 0.59 | 0.78 | 0.69 | 0.82 | 0.78 |
| StarbugII   | **0.038** | **0.046** | **0.043** | **0.027** | **0.034** |
| DOLPHOT     | **0.070** | **0.064** | **0.067** | **0.043** | **0.056** |

(σ = 1.4826 × MAD in log-flux ratio about the median; ~12 % flux scatter at 0.05 dex.)

See `figures/flux_scatter_norm_<tool>_<FILT>.png` for histograms; `figures/flux_magbin_<tool>_<FILT>.png` for the brightness-binned bias.

**Take-aways:**
- **StarbugII flux agreement with photutils is excellent — 4 % scatter** in matched bright stars across all five filters. This is the closest pairwise photometric match in the benchmark.
- **DOLPHOT flux scatter is ~12 %**, also good given that DOLPHOT runs on per-exposure cal frames while photutils runs on per-exposure crf frames and the band-merge happens differently. Most of the 12 % is probably real astrometric / PSF model difference between the two paths, not noise.
- **Crowdsource has 0.6-0.8 dex scatter (factor 4-7)** — clearly mis-photometering many matches. Combined with the low recovery fraction this means crowdsource is picking up junk and assigning huge fluxes to it. The cause is almost certainly that the runs were forced to `nskyx=nskyy=0` because crowdsource's sky-model code path crashes in `build_sparse_matrix` (line 501, `colnorm[startidx + syloc[i]]` is array-valued — see `benchmark/README.md`). Without a sky model, crowdsource interprets the structured PAH/dust background as point sources. Fixing that bug upstream would likely lift crowdsource from "broken" to "competitive".

## 5. Reported uncertainties (matched sources, median relative error)

| filter | photutils ref | crowdsource | StarbugII | DOLPHOT |
|--------|---------------|-------------|-----------|---------|
| F187N  | 0.41 | 0.15 | 0.085 | 0.137 |
| F210M  | 0.075 | 0.038 | 0.017 | 0.023 |
| F335M  | 0.075 | 0.032 | 0.013 | 0.015 |
| F470N  | 0.59 | 0.15 | 0.103 | 0.053 |
| F480M  | 0.095 | 0.054 | 0.017 | 0.017 |

**Take-aways:**
- **photutils-iter3 reports the largest uncertainties** by a wide margin — about 2× DOLPHOT in F480M and ~6× in F210M. The iter3 propagation pulls in per-exposure scatter as `flux_err_prop`, which captures real frame-to-frame variation that the other tools' single-fit-on-mosaic approach does not see. So photutils is more pessimistic and arguably more honest, though for an end-user the smaller DOLPHOT/StarbugII errors might be more useful for SNR-based selection cuts.
- **DOLPHOT and StarbugII agree to within ~10 % on relative errors** — they're computing the same statistical quantity (PSF-fit Cramer-Rao type uncertainty on a single image) on similar inputs.
- The narrow-band filters (F187N, F470N) show much larger uncertainties everywhere because the SNR is intrinsically lower; this is consistent across tools.

## 6. Effect of structured background emission

The Sickle field is dominated by spatially structured PAH / Brackett-α
nebular emission with brightness contrasts of order 10× across a single
detector. None of this background is saturated; the question is whether
each tool can model the smooth-on-arcsecond-scales nebular component
without confusing peaks in it for point sources, and whether the
reported uncertainties incorporate the local background level rather
than only the per-pixel readnoise + Poisson term.

- **photutils iter3** estimates the local background per fit using
  `LocalBackground(inner=2*FWHM, outer=4*FWHM)` inside `PSFPhotometry`
  and uses the JWST-pipeline ERR extension as the per-pixel weighting.
  The local-background subtraction is what keeps it from chasing PAH
  ridges; the consequence is that uncertainties grow in regions of
  steeper background gradient because `flux_err_prop` includes the
  variance contribution from the local-background fit.
- **DOLPHOT** estimates the sky per pixel via `calcsky` (annulus
  10–25 px, k-sigma clipped) before photometry, then re-estimates per
  source via `FitSky=2`. It does **not** read the JWST ERR extension —
  it builds its own per-pixel noise from `GAIN`, `EXPTIME`, and the
  read-noise number stored in the FITS header, plus a multiplicative
  buffer `NoiseMult=0.10`. So in regions where the JWST pipeline ERR
  is dominated by a flat-field / 1/f-correction component, DOLPHOT
  underestimates the noise and reports tighter error bars than warrant.
- **StarbugII** estimates the diffuse background via its own
  `-B` step (median-filter on a tunable kernel) and stores it in a
  `<base>-bgd.fits` map. The `-P` PSF photometry step does use the
  ERR extension. On Sickle this combination keeps PAH ridges out of
  the catalog and matches photutils on bright sources, but the median
  filter undercounts the smooth gradient between bright PAH peaks, so
  near-PAH faint stars are biased high.
- **Crowdsource** with the sky-parameter columns disabled (the original
  benchmark configuration, before the bug fix below) interpreted
  PAH-emission peaks as point sources — 25–35k tool-only "sources" per
  filter that have no photutils counterpart. With `nskyx=nskyy=0` the
  fitter has no smooth-component degree of freedom so it must put the
  background into stellar flux. After patching the bug at
  `crowdsource_base.py:501` (`syloc[i]` is array-valued; should be the
  scalar column index `i`) and re-running with `nskyx=nskyy=3`, the
  crowdsource recovery improves dramatically (see v2 results).

How this maps onto the metrics in §3-5:

- **Tool-only sources** (column `n_alt_only` in
  `extended_summary*.ecsv`) are dominated by PAH peaks for any tool
  that fails to fit a smooth background component.
- **Bias in flux ratios** is mostly local-background driven: tools
  that subtract a smaller local background produce higher fluxes and
  therefore a `median_flux_ratio` shifted from the photutils reference.
- **Reported uncertainty** is most honest in photutils because it
  includes per-frame propagation; DOLPHOT and StarbugII give
  Cramer-Rao bounds that ignore the variance contribution from the
  background-model uncertainty. With the iter3 pipeline computing
  `flux_err_prop` as the empirical frame-to-frame standard deviation,
  it captures the structured-background contribution that the other
  tools are missing.

Saturation is a separate axis of difficulty (driven by the bright OB
stars in the Quintuplet adjacent to the Sickle); each tool has its own
saturation handling, and the results are summarised in section 7.

## 7. Tuning / improvement notes per tool

### DOLPHOT
- **Astrometric registration in SW is the bottleneck** for recovery. Try (a) feeding DOLPHOT the iter3 catalog as a `xytfile` so it treats those as known stars rather than re-finding everything; (b) increasing `Align=4` -> `Align=6` or supplying explicit `img<i>_shift` from the JWST pipeline WCS.
- The default `RPSF=10` and `RAper=2.0` from the Weisz+24 NIRCam recipe are already in use; could tune `SecondPass=5` -> `7` for further faint-star recovery in dense regions.
- Residuals are clean: `dolphot/work/<FILT>/*.res.fits` shows no obvious stripes or ringing.
- Aperture correction is `ApCor=1` but reports `0 stars used` per image — DOLPHOT couldn't find clean isolated bright stars in these crowded subarrays. Run with an `ApCor=0` and apply zero-points externally, or pass an explicit aperture-correction file from a sparse-field calibration.

### crowdsource
- **Critical: patch `crowdsource_base.py:501`**. The bug `if colnorm[startidx + syloc[i]] == 0:` triggers when `nskyx, nskyy > 0` because `syloc[i]` is array-valued. Workaround `nskyx=nskyy=0` (used here) means no sky model and pollutes the catalog.
- After that fix, expect the recovery fraction to jump dramatically and the flux scatter to drop into DOLPHOT/StarbugII territory. Worth filing upstream.
- The Schlafly crowdsource pipeline has no JWST-specific PSF tools; we used the same stpsf grid PSF that the photutils pipeline does, wrapped in the same `WrappedPSFModel` class as in `crowdsource_catalogs_long.py`. So the PSF model is not the difference between crowdsource and the others.

### StarbugII
- Already produces good results out of the box. Two-step `-D -B -P` workflow with the auto-generated stpsf PSF library (4 × 4 = 16 NRCB1-4 PSFs + 1 NRCBLONG PSF per filter) was used.
- Three outputs per filter: `-ap` (aperture phot), `-psf` (PSF phot — used here), `-bgd` (background map).
- Could explore the `-A` artificial-star option for completeness curves.

### photutils iter3 (the reference)
- The iter3 cross-band union seeding is currently the most expensive part — the seed_union catalog has ~5k unique sources for Sickle. The per-exposure rerun with these seeds + `xy_bounds=±1 SW pix` (per memory `project_iter3_union_seeded.md`) fixed the systematic ~2× overfitting in F480M (memory `project_overfitting_root_cause.md`).
- Reported relative errors are larger than the other tools — these are propagated frame-to-frame, not single-image Cramer-Rao. That's a feature, but if external users compare error bars they should know this is comparing apples (per-frame propagated) to oranges (single-fit).
- The ~88k unique cross-band sources include narrow-band and saturation flags; the comparison here used `flux_<filt>` columns without applying the `near_saturated_<filt>` cut. Applying that cut would shift the matched-fraction numbers slightly upward.

## 8. Summary recommendation

For Sickle (and likely any other PAH-rich GC field):

| ranking | tool | comments |
|---------|------|----------|
| 1 | **photutils iter3** | most thorough, most sources, propagated errors, best satstar handling. |
| 2 | **StarbugII** | closest agreement with iter3 on matched sources (4 % flux), very good recovery. Best-in-class for "out-of-the-box JWST photometry". |
| 3 | **DOLPHOT** | excellent flux precision and residuals, but SW recovery limited by astrometric registration; conservative detection. |
| 4 | **crowdsource** | broken until the `build_sparse_matrix` shape bug is patched; with that fix, expected to be competitive with StarbugII. |

For an end user starting fresh on a NIRCam field with extended emission and **no in-house pipeline**, StarbugII gives the best bang/buck. For a project that already has a photutils pipeline, switching to iter3 + cross-band union seeding (as Sickle now does) outperforms all three alternative tools on raw source count, while DOLPHOT and StarbugII serve as useful independent cross-checks for the bright-end photometry.

## Reproducibility

All code, catalogs, and figures are under `/orange/adamginsburg/jwst/sickle/benchmark/`.
- `shared/config.py`, `shared/bandmerge.py`, `shared/compare_tools.py`, `shared/extended_compare.py`, `shared/load_dolphot.py`
- DOLPHOT param files, .phot, .res.fits per filter: `dolphot/work/<FILT>/`
- per-tool catalogs: `<tool>/catalogs/`
- band-merged FITS: `<tool>/bandmerged/<tool>_bandmerged_nrcb.fits`
- summary tables: `writeup/tables/{comparison_summary,extended_summary,detection_counts}.ecsv`
- figures: `writeup/figures/*.png`