Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

update COADD_FIBERSTATUS to bitwise OR when all inputs are bad #2093

Merged
merged 3 commits into from Aug 18, 2023

Conversation

sbailey
Copy link
Contributor

@sbailey sbailey commented Aug 11, 2023

This PR fixes the COADD_FIBERSTATUS bit propagation bug reported in #2032 .

For good targets, COADD_FIBERSTATUS is unchanged = the bitwise AND of input FIBERSTATUS.

For targets where all input spectra are bad, this PR changes COADD_FIBERSTATUS to the bitwise OR of input FIBERSTATUS values. This fixes the case in #2032 where the input spectra were all bad, but for different reasons, resulting in a bitwise OR of 0, incorrectly implying that the input FIBERSTATUS were ok. This will also result in other bad targets getting additional bits set (due to OR instead of AND).

Note that the incorrect COADD_FIBERSTATUS=0 cases had all ivar=0 anyway, which leads to ZWARN!=0 so we weren't getting incorrectly "good" redshifts for these anyway, but some analyses/bookkeeping were treating COADD_FIBERSTATUS==0 as "good hardware" so this makes it a bit more accurate.

import numpy as np
from desispec.io import read_spectra
import desispec.coaddition

sp = read_spectra('/global/cfs/cdirs/desi/spectro/redux/fuji/healpix/sv3/dark/281/28116/spectra-sv3-dark-28116.fits')
sp = sp[sp.fibermap['TARGETID'] == 39628362670410135]
print('Input spectra fibermap:')
print(sp.fibermap['TARGETID', 'FIBERSTATUS'])
print()

# make a copy since coaddition operates in-place
coadd = sp[:]
desispec.coaddition.coadd(coadd)
print('Coadded fibermap:')
print(coadd.fibermap['TARGETID', 'COADD_FIBERSTATUS'])

Results in

Input spectra fibermap:
     TARGETID     FIBERSTATUS
----------------- -----------
39628362670410135         256
39628362670410135           4

Coadded fibermap:
     TARGETID     COADD_FIBERSTATUS
----------------- -----------------
39628362670410135               260

(previously that case was incorrectly 0).

@ashleyjross @stephjuneau @akremin heads up. Comments welcome.

@stephjuneau
Copy link
Contributor

I haven't looked at the new code but from my recollection of the original code, the subtlety was that the fiberstatus needs to be determined on a per-amp basis in the case of the "badamp" bits and otherwise on a per spectrum basis. It might be worth to specifically find a case where only 1 or 2 amps are bad and the rest of the spectrum still useable for the unit tests.

@sbailey
Copy link
Contributor Author

sbailey commented Aug 18, 2023

Good question @stephjuneau . I reviewed N>>1 FIBERSTATUS cases with @akremin for what we really want COADD_FIBERSTATUS to mean and expanded the unit tests, which caught a problem in the original version of this PR for the case where 1 input spectrum is flagged bad for a per-fiber reason (like BROKENFIBER) while another input spectrum for the same target is flagged bad for a per-camera reason (like BADAMPB). Summarizing the final state:

If any input spectra are good:
    COADD_FIBERSTATUS = bitwise AND of the FIBERSTATUS of good input spectra, i.e.
                      = the AND of the FIBERSTATUSes actually used in the coadd
else (all are bad):
    COADD_FIBERSTATUS = bitwise OR of FIBERSTATUS input spectra

My first implementation added the else clause (bitwise OR), but was still taking the bitwise AND of all input spectra not just the ones actually used in the coadd.

Note: COADD_NUMEXP is the number of unique exposure IDs that contributed to the coadd, even if some cameras used fewer exposures than that due to BADAMPB/R/Z. I considered also implementing #2092 in this PR to add a record of which exposures were used for each camera, but I'm going to save that for a separate PR to avoid feature creep here. In any case, COADD_NUMEXP is only an approximately useful number since per-wavelength masking makes it a bin-by-bin value not just a camera-by-camera value, and the propagated IVAR is what you should really be using if the difference matters since coadding two high S/N spectra is very different from coadding one high S/N spectrum with one low S/N spectrum even if both are COADD_NUMEXP=2.

Along the way I also fixed #2088 reported by @araichoor where slicing spectra coadded in memory would crash due to coadd.scores being a dict instead of a Table.

@sbailey sbailey merged commit d405bf9 into main Aug 18, 2023
24 checks passed
@sbailey sbailey deleted the coadd_fiberstatus branch August 18, 2023 22:14
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

2 participants