I followed multiple ideas to improve the flux measurements, many of which were based in improving the detectorMaps:

• add DoubleDetectorMap
• deblend when centroiding
• flag lines with bad pixels when photometering
• review line list
• signal-to-noise cut in fitting detectorMap

In the end, I got down to about xSys=0.016 ySys=0.037 in fitting the detectorMap using sky lines (with readLineList.exclusionRadius=0.2):

reduceExposure.adjustDetectorMap INFO: Final fit: chi2=54099.876549 dof=23522 xRMS=0.024463 yRMS=0.037982 (0.003280 nm) from 11773/15601 lines
reduceExposure.adjustDetectorMap INFO: Fit quality from reserved lines: chi2=614201.080218 xRMS=0.024387 yRMS=0.064839 (0.005599 nm) from 1733 lines (10.0%)
reduceExposure.adjustDetectorMap INFO: Softening errors by x=0.016433, y=0.036511 pixels (0.003153 nm) to yield chi^2/dof=1
reduceExposure.adjustDetectorMap INFO: Softened fit: chi2=23387.268754 dof=23522 xRMS=0.025453 yRMS=0.042447 (0.003665 nm) from 11773 lines
reduceExposure.adjustDetectorMap INFO: Softened fit quality from reserved lines: chi2=244358.657181 xRMS=0.024418 yRMS=0.064045 (0.005530 nm) from 1733 lines

The difference between xSys and ySys likely indicates that there's more to be gained, but I'm not sure where yet. (xSys =~ ySys on arc lines, so the problem is intrinsic to the sky.)

With that detectorMap, I then get much better photometry, though still not ideal: chi^2 for the worst line fluxes is about 30 over five fibers, and that's when I include a 1% systematic floor (although before this work, the worst chi^2 was tens of thousands). It's possible that small inaccuracies in the PSF model are causing extra noise, and I think adding aperture corrections is the next step: if the fidelity of the PSF varies as a function of fiber or position within a fiber (which could be due to the PSF model failing to follow real variation, or the PSF model variation not being accurate), then we need aperture corrections.

import numpy as np
from pfs.drp.stella import ArcLineSet
lines = ArcLineSet.readFits("arcLines-046588-r1.fits")
select = np.isin(lines.fiberId, (255, 401, 464, 525, 587)) & np.isfinite(lines.intensity) & np.isfinite(lines.intensityErr) & (lines.status == 0) & ~lines.flag
intensityErr = np.where(lines.intensityErr > 0.01*lines.intensity, lines.intensityErr, 0.01*lines.intensity)
mean = {}
chi2 = {}
for wl in set(lines.wavelength[select]):
    choose = select & (lines.wavelength == wl)
    weight = 1.0/intensityErr[choose]**2
    mean[wl] = np.sum(lines.intensity[choose]*weight)/weight.sum()
    chi2[wl] = (((lines.intensity[choose] - mean[wl])/intensityErr[choose])**2).sum()

worst = sorted(chi2.keys(), key=lambda wl: chi2[wl], reverse=True)
print([(wl, chi2[wl]) for wl in worst[:10]])

yields the result:

[(753.27865, 30.999020699146175), (777.551, 30.257735938302954), (734.2907, 29.114637957334242), (885.22485, 25.031808119535704), (846.76845, 22.07350674889075), (691.45295, 21.767605415662246), (737.13955, 18.922802997897875), (854.10225, 18.811584139994014), (894.58505, 17.468980257836517), (961.0362, 16.9003032349872)]

Oh, I should mention that the OH doublets that we combined into single lines in our linelist are typically separated by something like a tenth (e.g., 6923.152+6923.287) to a few tenths (e.g., 7401.688+7402.029) of an Angstrom, which is about a tenth of a pixel to a few tenths of a pixel in arm=r, so maybe we shouldn't be surprised that there's an additional 0.03 pixel systematic error in the spectral dimension compared to the spatial. Centroiding on galaxies is noisier than centroiding on stars. Not sure if the measurement errors should account for that or not, but I don't think we need to worry too much.

Reviewed the changes in pfs_utils specifically to allow the repair of the integration test. Those specific changes are fine.

Now all changes approved. No major comments.

Merged.

I'm not sure that the blending is to blame. The raw centroid error will be a bit worse, but the code should be handling that (it's proportional to the FWHM/(SN)), but that should be very small for bright lines. The usual problem with galaxies is that the PSF changes, and that that changes the bias for the galaxy centroid. Is there any evidence of this?

Adding the deblender made very little difference, if any.

I'm not sure that I understand. You mean, simultaneous fitting of a double-line model with known splittings didn't improve the centroid repeatability?

No, I mean adding an SDSS-style deblender (removing pixels that don't belong to the peak of interest) didn't help.

I don't think I'd do that! This is much closer to forced-photometry part of the the crowded field problem, with low-spatial order tweaks to the "astrometry".