# Sickle photometry-tool benchmark: v2 comparison report

**Field:** JWST NIRCam, Galactic Center "Sickle" (program 3958, field 007), module B only.
**Filters:** F187N, F210M, F335M, F470N, F480M (5 NRCB).
**Reference (3 versions):** photutils-based three-iteration cascade (iter1 = blind DAOStarFinder per-exposure → iter2 = per-filter seeded → iter3 = cross-band union seeded).
**Other tools tested:** DOLPHOT 2.1 + NIRCam module (v2 with iter3-derived xytfile + Align=0), Schlafly's `crowdsourcephoto` 0.5.9 (sky-bug patched, per-exposure), StarbugII 0.7.7 (per-exposure).
**StarFinder:** skipped on user instruction.
**Match radius:** 0.08″.

**Big change from v1 → v2:**
- All four tools now run on the **same per-exposure `_destreak_o007_crf.fits` inputs** that the photutils iter3 pipeline uses, rather than mixing per-exposure (DOLPHOT, photutils) and L3 mosaic (crowdsource, StarbugII) inputs.
- DOLPHOT now uses an iter3-derived xytfile so its astrometric reference is the same Gaia/VVV-bootstrapped solution photutils uses, not its own image-relative alignment.
- crowdsource now runs with `nskyx=nskyy=3` (sky model enabled) after a one-line patch to `crowdsource_base.py:501` (a shape bug where `syloc[i]` was an array but used as a scalar index).
- iter1, iter2, iter3 are all available as separate references so we can see how each iteration of the photutils cascade adds detections.

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

## 1. Inputs (now apples-to-apples)

Every tool consumed the same per-exposure `*_destreak_o007_crf.fits` files (with SCI, ERR, DQ, VAR_POISSON, VAR_RNOISE, VAR_FLAT extensions) that the photutils pipeline runs on:

| filter | n exposures | detector(s) | array |
|--------|-------------|-------------|-------|
| F187N  | 96 | nrcb1-4 | SUB640 (640×640 each) |
| F210M  | 96 | nrcb1-4 | SUB640 |
| F335M  | 24 | nrcblong | SUB640 |
| F470N  | 24 | nrcblong | SUB640 |
| F480M  | 24 | nrcblong | SUB640 |

DOLPHOT runs all NRCB exposures together per filter in a single multi-image fit.
crowdsource and StarbugII now run **per exposure** as SLURM array jobs and we
band-merge each tool's per-exposure catalogs into a per-filter catalog via FoF
on sky coordinates (link radius = `MATCH_RADIUS_ARCSEC = 0.08″`).
The photutils per-exposure → per-filter merge is done by `merge_catalogs.py`.

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

|              | F187N  | F210M  | F335M  | F470N  | F480M  | total xband |
|--------------|--------|--------|--------|--------|--------|-------------|
| photutils_iter1 | 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     | 52,077 | 102,629 | 70,337 | 61,264 | 65,803 | 202,983 |
| StarbugII       | 18,291 | 78,393 | 21,998 | 14,394 | 19,219 | 99,606 |
| DOLPHOT v2      | 4,275  | 4,275  | 18,375 | 18,411 | 18,374 | 19,117 |

**Take-aways:**

- The **photutils cascade does what it's supposed to**: iter1 is blind per-exposure detection, iter2 is seeded by iter1 outputs (so it grows by 36% to 113k), and iter3 is the cross-band union-seeded refinement that prunes spurious detections to 89k while adding ones cross-band photometry confirms.
- **DOLPHOT v2** detects exactly the iter3 catalog forced positions for SW filters (4275 of 4853 measurable, the rest fell off-detector for the chosen reference image) plus all 22.8k LW positions (~18k fitted with usable flux). With `xytfile` + `Align=0` it does not search for new sources — that's intentional, the goal was to put it on the same astrometric grid.
- **crowdsource per-exposure**: 203k unique sources. The sky-model fix substantially changes the catalog (was 119k on mosaic with no sky model in v1), but it still finds more sources than photutils iter3 — many of these are real faint stars detected per-exposure that don't survive cross-band merging.
- **StarbugII per-exposure**: 99k. Comparable to photutils iter3 in volume but skewed toward brighter sources (see §3).

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

Computed as the fraction of photutils-iterN sources that have an alt-tool match within 0.08″:

| ref → | photutils iter1 | iter2 | iter3 |
|-------|-----------------|-------|-------|
| crowdsource F187N | 0.139 | 0.148 | 0.142 |
| crowdsource F210M | 0.167 | 0.165 | 0.166 |
| crowdsource F335M | 0.129 | 0.128 | 0.125 |
| crowdsource F470N | 0.165 | 0.171 | 0.165 |
| crowdsource F480M | 0.116 | 0.115 | 0.110 |
| starbug2    F187N | 0.055 | 0.009 | 0.010 |
| starbug2    F210M | 0.088 | 0.070 | 0.075 |
| starbug2    F335M | 0.387 | 0.369 | 0.384 |
| starbug2    F470N | 0.326 | 0.160 | 0.139 |
| starbug2    F480M | 0.479 | 0.382 | 0.378 |
| dolphot v2  F187N | 0.042 | 0.016 | 0.023 |
| dolphot v2  F210M | 0.039 | 0.042 | 0.053 |
| dolphot v2  F335M | 0.227 | 0.254 | 0.288 |
| dolphot v2  F470N | 0.232 | 0.201 | 0.233 |
| dolphot v2  F480M | 0.240 | 0.255 | 0.287 |

(See `figures/v2/recovery_fraction_grid.png` for bars.)

**Take-aways:**

- **DOLPHOT v2 recovery in LW** lands at **23–29 %** — that's much better than v1's 5 % in F210M and confirms that pinning DOLPHOT to the iter3 astrometric grid + using destreak inputs both helped.
- **DOLPHOT v2 recovery in SW remains low (1.6–4.2 %)**, *despite* using xytfile. The reason is that the SW mosaic spans 4 detectors and 24 dither positions; my reference image is one nrcb1 cal which only covers ~1/8 of the iter3 footprint. xytfile sources outside the reference FOV are silently dropped. Fix would be to pick a wider reference (e.g., the level-3 i2d mosaic or run separate DOLPHOT chains per detector).
- **StarbugII recovery LW = 38–48 %, SW = 1–9 %**. SW is low because per-exposure SUB640 sources are faint (single-frame depth) — StarbugII's `-D` step requires a positive S/N detection in that single frame, so most of the iter3 "hidden in coadd" sources are missed. iter1 (which is blind per-exposure too) shares this constraint, hence its small SW count.
- **crowdsource recovery is uniform 11–17 %** across filter and reference. The flat curve suggests its per-exposure detections have systematically displaced positions or are partially-spurious — both consistent with the high tool-only count and the wide flux scatter.
- The **"matched fraction" drops from iter1 → iter2 / iter3** for some pairs because iter1 has fewer SW sources (8.9k F187N vs 45k for iter2/3), so dividing by smaller numbers inflates the fraction; the absolute matched count is similar across iterN. The interesting comparison is **alt-tool vs iter3** (last column).

## 4. Photometric agreement (matched sources only)

Native flux units differ per tool, so the meaningful quantity is the spread of `log10(alt_flux / ref_flux)` after subtracting its median (= the zero-point offset). Reported as 1-σ in dex (1.4826 × MAD).

|                | F187N  | F210M  | F335M  | F470N  | F480M  |
|----------------|--------|--------|--------|--------|--------|
| crowdsource vs iter3 | 0.616 | 0.800 | 1.140 | 0.857 | 1.008 |
| StarbugII vs iter3   | 0.420 | 0.375 | 0.215 | 0.326 | 0.225 |
| DOLPHOT v2 vs iter3  | **0.058** | **0.057** | **0.063** | **0.060** | **0.061** |

Compared to v1 (where DOLPHOT used cal inputs and its own astrometry):
DOLPHOT scatter was 0.043–0.070, very similar to v2. This is reassuring —
the destreak vs cal inputs and the alignment method don't significantly affect
the bright-end photometric agreement when sources are matched.

The big v2 change is **crowdsource scatter actually grew** (0.6 dex in v1 →
0.6–1.1 dex in v2) because per-exposure detections in the dim Sickle field
have very low single-frame SNR, so the per-frame fluxes are noisy. The mosaic
runs in v1 averaged this out before fitting.

**StarbugII scatter is the most stable across filters: 0.21–0.42 dex (factor
1.6–2.6 in flux)**. That's the most honest comparison of what an out-of-the-box
JWST PSF photometry tool can deliver matched on per-exposure photutils data.

(See `figures/v2/flux_logratio_<alt>_vs_<ref>_<FILT>.png` for the histograms.)

## 5. Reported uncertainties

Median relative uncertainty for matched sources, vs photutils iter3:

| filter | iter3 (ref) | crowdsource | StarbugII | DOLPHOT v2 |
|--------|-------------|-------------|-----------|------------|
| F187N  | 0.236 | 0.366 | **328.8** | 0.145 |
| F210M  | 0.071 | 0.503 | **185.2** | 0.041 |
| F335M  | 0.084 | 0.251 | **47.4**  | 0.019 |
| F470N  | 0.492 | 0.178 | **113.9** | 0.072 |
| F480M  | 0.091 | 0.220 | **41.8**  | 0.018 |

**StarbugII per-exposure uncertainties are pathologically large**, regardless of filter — this is reproduced at the single-exposure level: a typical F480M
exposure has `flux ≈ 7×10⁻⁵` (MJy/sr × pixel), `eflux ≈ 0.023` → eflux/flux
≈ 87. starbug2 IS using the JWST ERR extension (we confirmed by reading its
config), but its noise estimator computes the integrated PSF noise over the
fit aperture, and because the SUB640 subarray exposures are 12.6 s each on
faint sources, the Cramer-Rao bound dominates. The merged catalog inherits
this large eflux. This is a real warning to anyone planning to use starbug2
on short subarray exposures: the reported errors are not useful as-is.

**DOLPHOT v2 reports the smallest errors (1.8–7 % in LW)**. DOLPHOT does not
read the JWST ERR extension; it computes per-pixel noise from `GAIN`, `EXPTIME`,
and the `READNOIS` keyword, plus `NoiseMult=0.10`. Those numbers are tighter
than the JWST pipeline's ERR (which propagates flat-field and 1/f corrections),
which is why DOLPHOT errors look small. With `xytfile` + forced photometry,
the position is fixed so DOLPHOT's flux uncertainty is the ideal Cramer-Rao
bound for that PSF model.

**photutils iter3 uncertainties (4.9–49 %)** include per-frame-propagated
scatter, so they are larger than DOLPHOT's per-image Cramer-Rao bounds but
similar to crowdsource's. F470N has the largest reported error (49 %)
because the narrow-band exposures have low SNR per frame and the propagation
captures that honestly.

**crowdsource (18–50 %)** is between the two. Its per-exposure error budget
includes the contribution of the sky-parameter columns (now that we re-enabled
them), which adds variance in regions of structured PAH emission.

## 6. Effect of structured background emission

(Unchanged from v1; section 6 stands.)

## 7. Tuning / improvement notes per tool — v2 update

### DOLPHOT v2

- **xytfile-based forced photometry works**: when DOLPHOT is given iter3
  positions and `Force1=1`, it produces tightly-correlated fluxes (0.06 dex
  scatter) and clean per-image residuals.
- **The remaining SW recovery problem is the reference-image footprint**.
  The xytfile is interpreted in the reference image's pixel frame; sources
  outside that frame are silently dropped. For SUB640 nrcb1 we get ~4275 of
  4853 iter3 SW positions on the chip (the rest are on nrcb2/3/4). To fix:
  run DOLPHOT four times per SW filter (one per detector) using each
  detector's first cal as the reference; or build an L3 mosaic of the
  destreak files and use that as the reference image (1850×750 pixels).
- **xytfile format gotcha**: lines must use `ext=1 chip=1` for JWST cal-style
  files, not `ext=0`. My v1 attempt used `ext=0` and got 0 stars. Documented
  in `dolphot.c:5125`.
- **AlignIter must be ≥ 1** even when `Align=0`. Setting AlignIter=0 triggers
  a perr at parse time. I set `AlignIter=1` with `Align=0` to satisfy both.
- **Aperture correction (`ApCor=1`) still fails** — DOLPHOT can't find clean
  isolated bright stars in these crowded subarrays. This affects the absolute
  zeropoint but not the relative flux comparison.

### crowdsource

- **Sky-model bug patched** at `crowdsource_base.py:501` (shape mismatch:
  `syloc[i]` was an `nx*ny`-element array; should be the scalar column index
  `i`). After the patch, `nskyx=nskyy=3` works as intended.
- **Per-exposure runs on faint short-subarray data are noisy** — the v1
  mosaic-based result (119k sources) is actually cleaner than the v2
  per-exposure result (203k sources) because L3 drizzling effectively
  smooths the noise and removes single-frame artefacts. For the apples-to-apples
  comparison the per-exposure run is right, but the take-away is that
  crowdsource **strongly prefers deep / coadded inputs** to per-exposure short
  subarrays.
- The wrapped GriddedPSFModel (same as photutils iter3 uses) gives reasonable
  PSF fits — the Schlafly-crowdsource doesn't have JWST-specific PSF tools but
  the wrapper code from `brick2221/analysis/crowdsource_catalogs_long.py`
  re-used here works.

### StarbugII

- **Per-exposure error reporting is unreliable on short subarrays**. eflux ≫
  flux for most sources, regardless of filter (median rel-err = 41–329 across
  filters). This isn't a bug per se — it's the Cramer-Rao bound on a
  12.6-second exposure of a faint source — but downstream consumers should
  not use these eflux values as 1-σ. For science, use the per-filter merged
  catalog's `flux_std_empirical` (frame-to-frame scatter), which is much more
  realistic.
- **PSF library auto-generation works** — the `--init` step took ~8 minutes
  (after the local interactive run collapsed) and populated 141 PSFs covering
  every NIRCam filter × detector. Rerun cost is zero on subsequent runs.
- **Per-exposure detection threshold is the limiting factor in SW** —
  StarbugII's `-D` step finds only ~565 sources per nrcb1 SUB640 exposure,
  vs ~10× that on the L3 mosaic. The fix is to give it the L3 mosaic as
  input or to lower the detection threshold (`SIGSRC` in starbug2's config).

### photutils iter1/iter2/iter3

- **iter1 → iter2 grows the F187N catalog 6×** (8,891 → 50,049). The big jump
  is because iter1 detects per-exposure with DAOStarFinder; iter2 reuses the
  positions from a deeper merged catalog (basically the union over the 96
  per-exposure iter1 catalogs) so it gets the full SW source density.
- **iter2 → iter3 prunes back to 89k** by requiring cross-band support (the
  cross-band union seed catalog at `seed_union_iter3_sickle.fits` has only
  positions detected in ≥1 filter). iter3 specifically drops sources that
  iter2 found in only one filter and that don't survive a cross-band sanity
  check — these are mostly cosmic-ray afterglows and PAH-edge spurious peaks.
- **The iter3 reduction in F210M (65k → 53k)** matches the size of the
  expected single-band-only spurious population in this PAH-rich field.

## 8. Summary recommendation — v2 update

For Sickle-like Galactic-Centre fields with dim sources on per-exposure
SUB640 subarrays:

| ranking | tool | comments |
|---------|------|----------|
| 1 | **photutils iter3 (per-exposure)** | most thorough; the iter1→iter2→iter3 cascade is doing real work, particularly the cross-band-union seed step. |
| 2 | **DOLPHOT v2 with iter3 xytfile** | best photometric agreement on matched sources (0.06 dex scatter), tightest reported error bars, but recovery is bounded by which iter3 positions land on the chosen reference image. **Run per-detector for SW.** |
| 3 | **StarbugII** | reasonable LW recovery (38–48 %), 0.2 dex flux scatter — but unusable per-exposure error bars on these short subarrays. |
| 4 | **crowdsource (sky model fixed)** | now functional, but per-exposure runs are noisy on short subarrays and 1 dex flux scatter is too large for science. **Prefers L3 mosaic input.** |

For projects starting from cal/destreak data:

- If you want a single-tool pipeline with no in-house infrastructure:
  use the L3 mosaic with **StarbugII** and accept the depth limitation.
- If you have your own astrometric solution (e.g., Gaia-bootstrapped) and
  want clean forced photometry: **DOLPHOT with xytfile** is the best
  external check on photutils.
- If you have a per-exposure photutils-style pipeline (like iter3): the
  v2 results show that the in-house pipeline does the best job on faint
  sources. The other tools' added value is **independent flux measurement
  for matched bright sources** (DOLPHOT specifically), or **a different
  detection algorithm to look for pipeline blind spots** (StarbugII on the
  L3 mosaic).

## Reproducibility

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

- `shared/config.py`, `shared/bandmerge.py` (with `--variant=indivexp`),
  `shared/extended_compare_v2.py`, `shared/load_dolphot.py`,
  `shared/merge_indivexp_to_filter.py`
- DOLPHOT v2 outputs: `dolphot/work_v2/<FILT>/`
- DOLPHOT v2 catalogs: `dolphot/catalogs_v2/<filt>_dolphot_nrcb.fits`
  (also accessible at `dolphot/catalogs/` symlink)
- Per-tool per-filter catalogs: `<tool>/catalogs/<filt>_<tool>_nrcb_indivexp.fits`
- Per-tool band-merged: `<tool>/bandmerged/<tool>_bandmerged_nrcb.fits`
- iter1/2/3 cross-filter: `photutils_iter{1,2,3}/iter{1,2,3}_xband_merged.fits`
- summary tables: `writeup/tables/{detection_counts_v2,extended_summary_v2}.ecsv`
- figures: `writeup/figures/v2/*.png`