Refactor local background and default chromatic corrections.

The local background subtraction code has been refactored to move the
computations into the fgcm code.  This allows fgcm to compute background offset
statistics and account for any net offsets in the zeropoints.  The chromatic
correction code has been updated to report the average chromatic correction per
ccd.  This yields better repeatability when the full pipeline does not include
chromatic terms (as is the case currently).  Additional tests have been added
to check for remaining net offsets and ensure they are small (<2 mmag).

python/lsst/fgcmcal/fgcmBuildStars.py

@@ -233,13 +233,12 @@ def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat,
         ccdKey = outputSchema['ccd'].asKey()
         instMagKey = outputSchema['instMag'].asKey()
         instMagErrKey = outputSchema['instMagErr'].asKey()
+        deltaMagBkgKey = outputSchema['deltaMagBkg'].asKey()
 
         # Prepare local background if desired
         if self.config.doSubtractLocalBackground:
             localBackgroundFluxKey = sourceSchema[self.config.localBackgroundFluxField].asKey()
             localBackgroundArea = np.pi*calibFluxApertureRadius**2.
-        else:
-            localBackground = 0.0
 
         aperOutputSchema = aperMapper.getOutputSchema()
 
@@ -271,11 +270,20 @@ def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat,
                 ccdId = dataRef.dataId[self.config.ccdDataRefName]
 
                 sources = dataRef.get(datasetType='src', flags=afwTable.SOURCE_IO_NO_FOOTPRINTS)
+                goodSrc = self.sourceSelector.selectSources(sources)
 
-                # If we are subtracting the local background, then correct here
-                # before we do the s/n selection.  This ensures we do not have
-                # bad stars after local background subtraction.
+                tempCat = afwTable.BaseCatalog(fullCatalog.schema)
+                tempCat.reserve(goodSrc.selected.sum())
+                tempCat.extend(sources[goodSrc.selected], mapper=sourceMapper)
+                tempCat[visitKey][:] = visit['visit']
+                tempCat[ccdKey][:] = ccdId
 
+                # Compute "instrumental magnitude" by scaling flux with exposure time.
+                scaledInstFlux = (sources[instFluxKey][goodSrc.selected] *
+                                  visit['scaling'][ccdMapping[ccdId]])
+                tempCat[instMagKey][:] = (-2.5*np.log10(scaledInstFlux) + 2.5*np.log10(expTime))
+
+                # Compute the change in magnitude from the background offset
                 if self.config.doSubtractLocalBackground:
                     # At the moment we only adjust the flux and not the flux
                     # error by the background because the error on
@@ -289,20 +297,14 @@ def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat,
                     # than a mmag in quadrature).
 
                     localBackground = localBackgroundArea*sources[localBackgroundFluxKey]
-                    sources[instFluxKey] -= localBackground
 
-                goodSrc = self.sourceSelector.selectSources(sources)
-
-                tempCat = afwTable.BaseCatalog(fullCatalog.schema)
-                tempCat.reserve(goodSrc.selected.sum())
-                tempCat.extend(sources[goodSrc.selected], mapper=sourceMapper)
-                tempCat[visitKey][:] = visit['visit']
-                tempCat[ccdKey][:] = ccdId
-
-                # Compute "instrumental magnitude" by scaling flux with exposure time.
-                scaledInstFlux = (sources[instFluxKey][goodSrc.selected] *
-                                  visit['scaling'][ccdMapping[ccdId]])
-                tempCat[instMagKey][:] = (-2.5*np.log10(scaledInstFlux) + 2.5*np.log10(expTime))
+                    # This is the difference between the mag with background correction
+                    # and the mag without background correction.
+                    tempCat[deltaMagBkgKey][:] = (-2.5*np.log10(sources[instFluxKey][goodSrc.selected] -
+                                                                localBackground[goodSrc.selected]) -
+                                                  -2.5*np.log10(sources[instFluxKey][goodSrc.selected]))
+                else:
+                    tempCat[deltaMagBkgKey][:] = 0.0
 
                 # Compute instMagErr from instFluxErr/instFlux, any scaling
                 # will cancel out.

python/lsst/fgcmcal/fgcmBuildStarsBase.py

@@ -599,6 +599,8 @@ def _makeSourceMapper(self, sourceSchema):
             "instMagErr", type=np.float32, doc="Instrumental magnitude error")
         sourceMapper.editOutputSchema().addField(
             "jacobian", type=np.float32, doc="Relative pixel scale from wcs jacobian")
+        sourceMapper.editOutputSchema().addField(
+            "deltaMagBkg", type=np.float32, doc="Change in magnitude due to local background offset")
 
         return sourceMapper
 

python/lsst/fgcmcal/fgcmBuildStarsTable.py

@@ -201,8 +201,6 @@ def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat,
         # Prepare local background if desired
         if self.config.doSubtractLocalBackground:
             localBackgroundArea = np.pi*calibFluxApertureRadius**2.
-        else:
-            localBackground = 0.0
 
         # Determine which columns we need from the sourceTable_visit catalogs
         columns = self._get_sourceTable_visit_columns()
@@ -221,10 +219,6 @@ def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat,
 
             df = srcTable.toDataFrame(columns)
 
-            if self.config.doSubtractLocalBackground:
-                localBackground = localBackgroundArea*df[self.config.localBackgroundFluxField].values
-                df[self.config.instFluxField] -= localBackground
-
             goodSrc = self.sourceSelector.selectSources(df)
             use, = np.where(goodSrc.selected)
 
@@ -239,6 +233,28 @@ def fgcmMakeAllStarObservations(self, groupedDataRefs, visitCat,
             tempCat[ccdKey][:] = df[self.config.ccdDataRefName].values[use]
             tempCat['psf_candidate'] = df['Calib_psf_candidate'].values[use]
 
+            if self.config.doSubtractLocalBackground:
+                # At the moment we only adjust the flux and not the flux
+                # error by the background because the error on
+                # base_LocalBackground_instFlux is the rms error in the
+                # background annulus, not the error on the mean in the
+                # background estimate (which is much smaller, by sqrt(n)
+                # pixels used to estimate the background, which we do not
+                # have access to in this task).  In the default settings,
+                # the annulus is sufficiently large such that these
+                # additional errors are are negligibly small (much less
+                # than a mmag in quadrature).
+
+                localBackground = localBackgroundArea*df[self.config.localBackgroundFluxField].values
+
+                # This is the difference between the mag with local background correction
+                # and the mag without local background correction.
+                tempCat['deltaMagBkg'] = (-2.5*np.log10(df[self.config.instFluxField].values[use] -
+                                                        localBackground[use]) -
+                                          -2.5*np.log10(df[self.config.instFluxField].values[use]))
+            else:
+                tempCat['deltaMagBkg'][:] = 0.0
+
             # Need to loop over ccds here
             for detector in camera:
                 ccdId = detector.getId()

python/lsst/fgcmcal/fgcmCalibrateTractBase.py

@@ -325,6 +325,7 @@ def runDataRef(self, butler, dataRefs):
                             fgcmStarIdCat['nObs'][:],
                             obsX=fgcmStarObservationCat['x'][obsIndex],
                             obsY=fgcmStarObservationCat['y'][obsIndex],
+                            obsDeltaMagBkg=fgcmStarObservationCat['deltaMagBkg'][obsIndex],
                             psfCandidate=fgcmStarObservationCat['psf_candidate'][obsIndex],
                             refID=refId,
                             refMag=refMag,

python/lsst/fgcmcal/fgcmFitCycle.py

@@ -231,6 +231,18 @@ class FgcmFitCycleConfig(pexConfig.Config):
         dtype=int,
         default=50,
     )
+    deltaMagBkgOffsetPercentile = pexConfig.Field(
+        doc=("Percentile brightest stars on a visit/ccd to use to compute net "
+             "offset from local background subtraction."),
+        dtype=float,
+        default=0.25,
+    )
+    deltaMagBkgPerCcd = pexConfig.Field(
+        doc=("Compute net offset from local background subtraction per-ccd? "
+             "Otherwise, use computation per visit."),
+        dtype=bool,
+        default=False,
+    )
     utBoundary = pexConfig.Field(
         doc="Boundary (in UTC) from day-to-day",
         dtype=float,
@@ -920,6 +932,7 @@ def _fgcmFitCycle(self, butler):
                             starIds['nObs'][:],
                             obsX=starObs['x'][starIndices['obsIndex']],
                             obsY=starObs['y'][starIndices['obsIndex']],
+                            obsDeltaMagBkg=starObs['deltaMagBkg'][starIndices['obsIndex']],
                             psfCandidate=starObs['psf_candidate'][starIndices['obsIndex']],
                             refID=refId,
                             refMag=refMag,
@@ -1111,6 +1124,8 @@ def _loadParameters(self, butler):
                                        (parCat['compMirrorChromaticity'].size, )),
                                       ('MIRRORCHROMATICITYPIVOT', 'f8',
                                        (parCat['mirrorChromaticityPivot'].size, )),
+                                      ('COMPMEDIANSEDSLOPE', 'f8',
+                                       (parCat['compMedianSedSlope'].size, )),
                                       ('COMPAPERCORRPIVOT', 'f8',
                                        (parCat['compAperCorrPivot'].size, )),
                                       ('COMPAPERCORRSLOPE', 'f8',
@@ -1131,6 +1146,8 @@ def _loadParameters(self, butler):
                                        (parCat['compExpGray'].size, )),
                                       ('COMPVARGRAY', 'f8',
                                        (parCat['compVarGray'].size, )),
+                                      ('COMPEXPDELTAMAGBKG', 'f8',
+                                       (parCat['compExpDeltaMagBkg'].size, )),
                                       ('COMPNGOODSTARPEREXP', 'i4',
                                        (parCat['compNGoodStarPerExp'].size, )),
                                       ('COMPSIGFGCM', 'f8',
@@ -1165,6 +1182,7 @@ def _loadParameters(self, butler):
         inParams['COMPREFSIGMA'][:] = parCat['compRefSigma'][0, :]
         inParams['COMPMIRRORCHROMATICITY'][:] = parCat['compMirrorChromaticity'][0, :]
         inParams['MIRRORCHROMATICITYPIVOT'][:] = parCat['mirrorChromaticityPivot'][0, :]
+        inParams['COMPMEDIANSEDSLOPE'][:] = parCat['compMedianSedSlope'][0, :]
         inParams['COMPAPERCORRPIVOT'][:] = parCat['compAperCorrPivot'][0, :]
         inParams['COMPAPERCORRSLOPE'][:] = parCat['compAperCorrSlope'][0, :]
         inParams['COMPAPERCORRSLOPEERR'][:] = parCat['compAperCorrSlopeErr'][0, :]
@@ -1175,6 +1193,7 @@ def _loadParameters(self, butler):
         inParams['COMPMODELERRPARS'][:] = parCat['compModelErrPars'][0, :]
         inParams['COMPEXPGRAY'][:] = parCat['compExpGray'][0, :]
         inParams['COMPVARGRAY'][:] = parCat['compVarGray'][0, :]
+        inParams['COMPEXPDELTAMAGBKG'][:] = parCat['compExpDeltaMagBkg'][0, :]
         inParams['COMPNGOODSTARPEREXP'][:] = parCat['compNGoodStarPerExp'][0, :]
         inParams['COMPSIGFGCM'][:] = parCat['compSigFgcm'][0, :]
         inParams['COMPSIGMACAL'][:] = parCat['compSigmaCal'][0, :]
@@ -1347,6 +1366,9 @@ def _makeParSchema(self, parInfo, pars, parSuperStarFlat,
         parSchema.addField('mirrorChromaticityPivot', type='ArrayD',
                            doc='Mirror chromaticity pivot mjd',
                            size=pars['MIRRORCHROMATICITYPIVOT'].size)
+        parSchema.addField('compMedianSedSlope', type='ArrayD',
+                           doc='Computed median SED slope (per band)',
+                           size=pars['COMPMEDIANSEDSLOPE'].size)
         parSchema.addField('compAperCorrPivot', type='ArrayD', doc='Aperture correction pivot',
                            size=pars['COMPAPERCORRPIVOT'].size)
         parSchema.addField('compAperCorrSlope', type='ArrayD', doc='Aperture correction slope',
@@ -1367,6 +1389,9 @@ def _makeParSchema(self, parInfo, pars, parSuperStarFlat,
                            size=pars['COMPEXPGRAY'].size)
         parSchema.addField('compVarGray', type='ArrayD', doc='Computed exposure variance',
                            size=pars['COMPVARGRAY'].size)
+        parSchema.addField('compExpDeltaMagBkg', type='ArrayD',
+                           doc='Computed exposure offset due to background',
+                           size=pars['COMPEXPDELTAMAGBKG'].size)
         parSchema.addField('compNGoodStarPerExp', type='ArrayI',
                            doc='Computed number of good stars per exposure',
                            size=pars['COMPNGOODSTARPEREXP'].size)
@@ -1444,7 +1469,7 @@ def _makeParCatalog(self, parSchema, parInfo, pars, parSuperStarFlat,
                     'compModelErrExptimePivot', 'compModelErrFwhmPivot',
                     'compModelErrSkyPivot', 'compModelErrPars',
                     'compExpGray', 'compVarGray', 'compNGoodStarPerExp', 'compSigFgcm',
-                    'compSigmaCal',
+                    'compSigmaCal', 'compExpDeltaMagBkg', 'compMedianSedSlope',
                     'compRetrievedLnPwv', 'compRetrievedLnPwvRaw', 'compRetrievedLnPwvFlag',
                     'compRetrievedTauNight']
 

python/lsst/fgcmcal/fgcmOutputProducts.py

@@ -99,6 +99,11 @@ class FgcmOutputProductsConfig(pexConfig.Config):
         dtype=bool,
         default=True,
     )
+    doApplyMeanChromaticCorrection = pexConfig.Field(
+        doc="Apply the mean chromatic correction to the zeropoints?",
+        dtype=bool,
+        default=True,
+    )
     refObjLoader = pexConfig.ConfigurableField(
         target=LoadIndexedReferenceObjectsTask,
         doc="reference object loader for 'absolute' photometric calibration",
@@ -834,12 +839,21 @@ def _outputZeropoints(self, butler, zptCat, visitCat, offsets, tract=None):
             # Retrieve overall scaling
             scaling = scalingMapping[rec['visit']][ccdMapping[rec['ccd']]]
 
+            # The postCalibrationOffset describe any zeropoint offsets
+            # to apply after the fgcm calibration.  The first part comes
+            # from the reference catalog match (used in testing).  The
+            # second part comes from the mean chromatic correction
+            # (if configured).
+            postCalibrationOffset = offsetMapping[rec['filtername']]
+            if self.config.doApplyMeanChromaticCorrection:
+                postCalibrationOffset += rec['fgcmDeltaChrom']
+
             fgcmSuperStarField = self._getChebyshevBoundedField(rec['fgcmfZptSstarCheb'],
                                                                 rec['fgcmfZptChebXyMax'])
             # Convert from FGCM AB to nJy
             fgcmZptField = self._getChebyshevBoundedField((rec['fgcmfZptCheb']*units.AB).to_value(units.nJy),
                                                           rec['fgcmfZptChebXyMax'],
-                                                          offset=offsetMapping[rec['filtername']],
+                                                          offset=postCalibrationOffset,
                                                           scaling=scaling)
 
             if self.config.doComposeWcsJacobian:

python/lsst/fgcmcal/utilities.py

@@ -126,6 +126,8 @@ def makeConfigDict(config, log, camera, maxIter,
                   'ccdGraySubCCDTriangular': config.ccdGraySubCcdTriangular,
                   'cycleNumber': config.cycleNumber,
                   'maxIter': maxIter,
+                  'deltaMagBkgOffsetPercentile': config.deltaMagBkgOffsetPercentile,
+                  'deltaMagBkgPerCcd': config.deltaMagBkgPerCcd,
                   'UTBoundary': config.utBoundary,
                   'washMJDs': config.washMjds,
                   'epochMJDs': config.epochMjds,
@@ -560,7 +562,10 @@ def makeZptSchema(superStarChebyshevSize, zptChebyshevSize):
                        doc='Retrieved i10 integral, estimated from stars (only for flag 1)')
     zptSchema.addField('fgcmGry', type=np.float64,
                        doc='Estimated gray extinction relative to atmospheric solution; '
-                       'only for flag <= 4')
+                       'only for fgcmFlag <= 4 (see fgcmFlag) ')
+    zptSchema.addField('fgcmDeltaChrom', type=np.float64,
+                       doc='Mean chromatic correction for stars in this ccd; '
+                       'only for fgcmFlag <= 4 (see fgcmFlag)')
     zptSchema.addField('fgcmZptVar', type=np.float64, doc='Variance of zeropoint over ccd')
     zptSchema.addField('fgcmTilings', type=np.float64,
                        doc='Number of photometric tilings used for solution for ccd')
@@ -587,6 +592,10 @@ def makeZptSchema(superStarChebyshevSize, zptChebyshevSize):
                        'at the time of the exposure')
     zptSchema.addField('fgcmFlat', type=np.float64, doc='Superstarflat illumination correction')
     zptSchema.addField('fgcmAperCorr', type=np.float64, doc='Aperture correction estimated by fgcm')
+    zptSchema.addField('fgcmDeltaMagBkg', type=np.float64,
+                       doc=('Local background correction from brightest percentile '
+                            '(value set by deltaMagBkgOffsetPercentile) calibration '
+                            'stars.'))
     zptSchema.addField('exptime', type=np.float32, doc='Exposure time')
     zptSchema.addField('filtername', type=str, size=10, doc='Filter name')
 
@@ -630,6 +639,7 @@ def makeZptCat(zptSchema, zpStruct):
     zptCat['fgcmR0'][:] = zpStruct['FGCM_R0']
     zptCat['fgcmR10'][:] = zpStruct['FGCM_R10']
     zptCat['fgcmGry'][:] = zpStruct['FGCM_GRY']
+    zptCat['fgcmDeltaChrom'][:] = zpStruct['FGCM_DELTACHROM']
     zptCat['fgcmZptVar'][:] = zpStruct['FGCM_ZPTVAR']
     zptCat['fgcmTilings'][:] = zpStruct['FGCM_TILINGS']
     zptCat['fgcmFpGry'][:] = zpStruct['FGCM_FPGRY']
@@ -641,6 +651,7 @@ def makeZptCat(zptSchema, zpStruct):
     zptCat['fgcmDust'][:] = zpStruct['FGCM_DUST']
     zptCat['fgcmFlat'][:] = zpStruct['FGCM_FLAT']
     zptCat['fgcmAperCorr'][:] = zpStruct['FGCM_APERCORR']
+    zptCat['fgcmDeltaMagBkg'][:] = zpStruct['FGCM_DELTAMAGBKG']
     zptCat['exptime'][:] = zpStruct['EXPTIME']
 
     return zptCat