Merge pull request #42 from lsst/tickets/DM-26085

DM-26085: Fix fgcmcal zeropoint offset due to background offset.

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