Merge pull request #84 from lsst/tickets/DM-31880

DM-31880: Add support for new fgcm reference star options.

python/lsst/fgcmcal/fgcmFitCycle.py

@@ -329,10 +329,17 @@ class FgcmFitCycleConfig(pipeBase.PipelineTaskConfig,
         default=4.0,
     )
     applyRefStarColorCuts = pexConfig.Field(
-        doc="Apply color cuts to reference stars?",
+        doc=("Apply color cuts defined in ``starColorCuts`` to reference stars? "
+             "These cuts are in addition to any cuts defined in ``refStarColorCuts``"),
         dtype=bool,
         default=True,
     )
+    useExposureReferenceOffset = pexConfig.Field(
+        doc=("Use per-exposure (visit) offsets between calibrated stars and reference stars "
+             "for final zeropoints? This may help uniformity for disjoint surveys."),
+        dtype=bool,
+        default=False,
+    )
     nCore = pexConfig.Field(
         doc="Number of cores to use",
         dtype=int,
@@ -659,7 +666,18 @@ class FgcmFitCycleConfig(pipeBase.PipelineTaskConfig,
         default=None,
     )
     starColorCuts = pexConfig.ListField(
-        doc="Encoded star-color cuts (to be cleaned up)",
+        doc=("Encoded star-color cuts (using calibration star colors). "
+             "This is a list with each entry a string of the format "
+             "``band1,band2,low,high`` such that only stars of color "
+             "low < band1 - band2 < high will be used for calibration."),
+        dtype=str,
+        default=("NO_DATA",),
+    )
+    refStarColorCuts = pexConfig.ListField(
+        doc=("Encoded star color cuts specifically to apply to reference stars. "
+             "This is a list with each entry a string of the format "
+             "``band1,band2,low,high`` such that only stars of color "
+             "low < band1 - band2 < high will be used as reference stars."),
         dtype=str,
         default=("NO_DATA",),
     )
@@ -1295,6 +1313,8 @@ def _loadParameters(self, parCat):
                                        (parCat['compExpDeltaMagBkg'].size, )),
                                       ('COMPNGOODSTARPEREXP', 'i4',
                                        (parCat['compNGoodStarPerExp'].size, )),
+                                      ('COMPEXPREFOFFSET', 'f8',
+                                       (parCat['compExpRefOffset'].size, )),
                                       ('COMPSIGFGCM', 'f8',
                                        (parCat['compSigFgcm'].size, )),
                                       ('COMPSIGMACAL', 'f8',
@@ -1354,6 +1374,7 @@ def _loadParameters(self, parCat):
         inParams['COMPVARGRAY'][:] = parCat['compVarGray'][0, :]
         inParams['COMPEXPDELTAMAGBKG'][:] = parCat['compExpDeltaMagBkg'][0, :]
         inParams['COMPNGOODSTARPEREXP'][:] = parCat['compNGoodStarPerExp'][0, :]
+        inParams['COMPEXPREFOFFSET'][:] = parCat['compExpRefOffset'][0, :]
         inParams['COMPSIGFGCM'][:] = parCat['compSigFgcm'][0, :]
         inParams['COMPSIGMACAL'][:] = parCat['compSigmaCal'][0, :]
         inParams['COMPRETRIEVEDLNPWV'][:] = parCat['compRetrievedLnPwv'][0, :]
@@ -1563,6 +1584,9 @@ def _makeParSchema(self, parInfo, pars, parSuperStarFlat,
         parSchema.addField('compNGoodStarPerExp', type='ArrayI',
                            doc='Computed number of good stars per exposure',
                            size=pars['COMPNGOODSTARPEREXP'].size)
+        parSchema.addField('compExpRefOffset', type='ArrayD',
+                           doc='Computed per-visit median offset between standard stars and ref stars.',
+                           size=pars['COMPEXPREFOFFSET'].size)
         parSchema.addField('compSigFgcm', type='ArrayD', doc='Computed sigma_fgcm (intrinsic repeatability)',
                            size=pars['COMPSIGFGCM'].size)
         parSchema.addField('compSigmaCal', type='ArrayD', doc='Computed sigma_cal (systematic error floor)',
@@ -1662,7 +1686,7 @@ def _makeParCatalog(self, parSchema, parInfo, pars, parSuperStarFlat,
                     'compRetrievedLnPwv', 'compRetrievedLnPwvRaw', 'compRetrievedLnPwvFlag',
                     'compRetrievedTauNight', 'compEpsilon', 'compMedDeltaAper',
                     'compGlobalEpsilon', 'compEpsilonMap', 'compEpsilonNStarMap',
-                    'compEpsilonCcdMap', 'compEpsilonCcdNStarMap']
+                    'compEpsilonCcdMap', 'compEpsilonCcdNStarMap', 'compExpRefOffset']
 
         for scalarName in scalarNames:
             rec[scalarName] = pars[scalarName.upper()]

python/lsst/fgcmcal/utilities.py

@@ -83,9 +83,21 @@ def makeConfigDict(config, log, camera, maxIter,
     # process the starColorCuts
     starColorCutList = []
     for ccut in config.starColorCuts:
+        if ccut == 'NO_DATA':
+            # No color cuts to apply.
+            break
         parts = ccut.split(',')
         starColorCutList.append([parts[0], parts[1], float(parts[2]), float(parts[3])])
 
+    # process the refStarColorCuts
+    refStarColorCutList = []
+    for ccut in config.refStarColorCuts:
+        if ccut == 'NO_DATA':
+            # No color cuts to apply.
+            break
+        parts = ccut.split(',')
+        refStarColorCutList.append([parts[0], parts[1], float(parts[2]), float(parts[3])])
+
     # TODO: Having direct access to the mirror area from the camera would be
     #  useful.  See DM-16489.
     # Mirror area in cm**2
@@ -169,8 +181,10 @@ def makeConfigDict(config, log, camera, maxIter,
                   'refStarSnMin': config.refStarSnMin,
                   'refStarOutlierNSig': config.refStarOutlierNSig,
                   'applyRefStarColorCuts': config.applyRefStarColorCuts,
+                  'useExposureReferenceOffset': config.useExposureReferenceOffset,
                   'illegalValue': FGCM_ILLEGAL_VALUE,  # internally used by fgcm.
                   'starColorCuts': starColorCutList,
+                  'refStarColorCuts': refStarColorCutList,
                   'aperCorrFitNBins': config.aperCorrFitNBins,
                   'aperCorrInputSlopeDict': dict(config.aperCorrInputSlopeDict),
                   'sedBoundaryTermDict': config.sedboundaryterms.toDict()['data'],

tests/config/fgcmFitCycleHsc.py

@@ -51,6 +51,8 @@
                                             extrapolated=False, constant=0.75)
 # Define good stars with an r-i color cut.
 config.starColorCuts = ('r,i,-0.50,2.25',)
+config.refStarColorCuts = ('r, i, 0.0, 2.0',)
+config.useExposureReferenceOffset = True
 config.precomputeSuperStarInitialCycle = False
 config.superStarSubCcdDict = {'g': True,
                               'r': True,

tests/test_fgcmcal_hsc.py

@@ -92,25 +92,27 @@ def test_fgcmcalPipeline(self):
         nOkZp = 27
         nBadZp = 1093
         nStdStars = 235
-        nPlots = 43
+        nPlots = 48
 
         self._testFgcmFitCycle(instName, testName,
                                0, nZp, nGoodZp, nOkZp, nBadZp, nStdStars, nPlots, skipChecks=True)
         self._testFgcmFitCycle(instName, testName,
                                1, nZp, nGoodZp, nOkZp, nBadZp, nStdStars, nPlots, skipChecks=True)
 
         # We need to create an extra config file to turn on "sub-ccd gray" for testing.
+        # We also want to exercise the code path setting useExposureReferenceOffset = False.
         extraConfigFile = os.path.join(self.testDir, "cycle03_patch_config.py")
         with open(extraConfigFile, "w") as f:
             f.write("config.isFinalCycle = True\n")
             f.write("config.ccdGraySubCcdDict = {'g': True, 'r': True, 'i': True}\n")
+            f.write("config.useExposureReferenceOffset = False")
 
         self._testFgcmFitCycle(instName, testName,
                                2, nZp, nGoodZp, nOkZp, nBadZp, nStdStars, nPlots,
                                extraConfig=extraConfigFile)
 
-        zpOffsets = np.array([-0.0008161436999216676,
-                              0.006149172317236662])
+        zpOffsets = np.array([-0.0008051003096625209,
+                              0.0072303167544305325])
 
         self._testFgcmOutputProducts(instName, testName,
                                      zpOffsets, 36236, 87, 'i', 1)
@@ -135,8 +137,8 @@ def test_fgcmcalMultipleFitPipeline(self):
 
         # These are slightly different from above due to the configuration change
         # mid-way in the separate fits.
-        zpOffsets = np.array([-0.0007102462113834918,
-                              0.005907602142542601])
+        zpOffsets = np.array([-0.0006988655077293515,
+                              0.004102597013115883])
 
         self._testFgcmMultiFit(instName, testName,
                                "physical_filter IN ('HSC-G', 'HSC-R', 'HSC-I')",
@@ -159,8 +161,8 @@ def test_fgcmcalTractPipeline(self):
                               nBand, i0Std, i0Recon, i10Std, i10Recon)
 
         rawRepeatability = np.array([0.0,
-                                     0.004436014222072738,
-                                     0.00451764656339253])
+                                     0.0025195920941720683,
+                                     0.004095912225403857])
         filterNCalibMap = {'HSC-R': 12,
                            'HSC-I': 15}