Merge pull request #96 from lsst/tickets/DM-34061

DM-34061: Add new outputs for standard atmosphere and passbands.

python/lsst/fgcmcal/fgcmMakeLut.py

@@ -39,6 +39,7 @@
 from lsst.pipe.base import connectionTypes
 import lsst.afw.table as afwTable
 import lsst.afw.cameraGeom as afwCameraGeom
+from lsst.afw.image import TransmissionCurve
 from .utilities import lookupStaticCalibrations
 
 import fgcm
@@ -98,6 +99,21 @@ class FgcmMakeLutConnections(pipeBase.PipelineTaskConnections,
         dimensions=("instrument",),
     )
 
+    fgcmStandardAtmosphere = connectionTypes.Output(
+        doc="Standard atmosphere used for FGCM calibration.",
+        name="fgcm_standard_atmosphere",
+        storageClass="TransmissionCurve",
+        dimensions=("instrument",),
+    )
+
+    fgcmStandardPassbands = connectionTypes.Output(
+        doc="Standard passbands used for FGCM calibration.",
+        name="fgcm_standard_passband",
+        storageClass="TransmissionCurve",
+        dimensions=("instrument", "physical_filter"),
+        multiple=True,
+    )
+
 
 class FgcmMakeLutParametersConfig(pexConfig.Config):
     """Config for parameters if atmosphereTableName not available"""
@@ -304,11 +320,18 @@ def runQuantum(self, butlerQC, inputRefs, outputRefs):
         filterHandleDict = {filterHandle.dataId['physical_filter']: filterHandle for
                             filterHandle in filterHandles}
 
-        lutCat = self._fgcmMakeLut(camera,
+        struct = self._fgcmMakeLut(camera,
                                    opticsHandle,
                                    sensorHandleDict,
                                    filterHandleDict)
-        butlerQC.put(lutCat, outputRefs.fgcmLookUpTable)
+
+        butlerQC.put(struct.fgcmLookUpTable, outputRefs.fgcmLookUpTable)
+        butlerQC.put(struct.fgcmStandardAtmosphere, outputRefs.fgcmStandardAtmosphere)
+
+        refDict = {passbandRef.dataId['physical_filter']: passbandRef for
+                   passbandRef in outputRefs.fgcmStandardPassbands}
+        for physical_filter, passband in struct.fgcmStandardPassbands.items():
+            butlerQC.put(passband, refDict[physical_filter])
 
     def _fgcmMakeLut(self, camera, opticsHandle, sensorHandleDict,
                      filterHandleDict):
@@ -330,8 +353,16 @@ def _fgcmMakeLut(self, camera, opticsHandle, sensorHandleDict,
 
         Returns
         -------
-        fgcmLookUpTable : `BaseCatalog`
-            The FGCM look-up table.
+        retStruct : `lsst.pipe.base.Struct`
+            Output structure with keys:
+
+            fgcmLookUpTable : `BaseCatalog`
+                The FGCM look-up table.
+            fgcmStandardAtmosphere : `lsst.afw.image.TransmissionCurve`
+                Transmission curve for the FGCM standard atmosphere.
+            fgcmStandardPassbands : `dict` [`str`, `lsst.afw.image.TransmissionCurve`]
+                Dictionary of standard passbands, with the key as the
+                physical filter name.
         """
         # number of ccds from the length of the camera iterator
         nCcd = len(camera)
@@ -393,7 +424,31 @@ def _fgcmMakeLut(self, camera, opticsHandle, sensorHandleDict,
 
         lutCat = self._makeLutCat(lutSchema, physicalFilterString,
                                   stdPhysicalFilterString, atmosphereTableName)
-        return lutCat
+
+        atmStd = TransmissionCurve.makeSpatiallyConstant(
+            throughput=self.fgcmLutMaker.atmStdTrans.astype(np.float64),
+            wavelengths=self.fgcmLutMaker.atmLambda.astype(np.float64),
+            throughputAtMin=self.fgcmLutMaker.atmStdTrans[0],
+            throughputAtMax=self.fgcmLutMaker.atmStdTrans[1],
+        )
+
+        fgcmStandardPassbands = {}
+        for i, physical_filter in enumerate(self.fgcmLutMaker.filterNames):
+            passband = self.fgcmLutMaker.throughputs[i]['THROUGHPUT_AVG']*self.fgcmLutMaker.atmStdTrans
+            fgcmStandardPassbands[physical_filter] = TransmissionCurve.makeSpatiallyConstant(
+                throughput=passband.astype(np.float64),
+                wavelengths=self.fgcmLutMaker.atmLambda.astype(np.float64),
+                throughputAtMin=passband[0],
+                throughputAtMax=passband[-1],
+            )
+
+        retStruct = pipeBase.Struct(
+            fgcmLookUpTable=lutCat,
+            fgcmStandardAtmosphere=atmStd,
+            fgcmStandardPassbands=fgcmStandardPassbands,
+        )
+
+        return retStruct
 
     def _getStdPhysicalFilterList(self):
         """Get the standard physical filter lists from config.physicalFilters
@@ -623,7 +678,7 @@ def _makeLutCat(self, lutSchema, physicalFilterString, stdPhysicalFilterString,
         Returns
         -------
         lutCat: `afwTable.BaseCatalog`
-           Lut catalog for persistence
+           Look-up table catalog for persistence.
         """
 
         # The somewhat strange format is to make sure that

tests/fgcmcalTestBase.py

@@ -209,6 +209,25 @@ def _testFgcmMakeLut(self, instName, testName, nBand, i0Std, i0Recon, i10Std, i1
 
         self.assertFloatsAlmostEqual(i10Recon, i1/i0, msg='i10Recon', rtol=1e-5)
 
+        # Check that the standard atmosphere was output and non-zero.
+        atmStd = butler.get('fgcm_standard_atmosphere',
+                            collections=[outputCollection],
+                            instrument=instName)
+        bounds = atmStd.getWavelengthBounds()
+        lambdas = np.linspace(bounds[0], bounds[1], 1000)
+        tputs = atmStd.sampleAt(position=geom.Point2D(0.0, 0.0), wavelengths=lambdas)
+        self.assertGreater(np.min(tputs), 0.0)
+
+        # Check that the standard passbands were output and non-zero.
+        for physical_filter in fgcmLut.filterNames:
+            passband = butler.get('fgcm_standard_passband',
+                                  collections=[outputCollection],
+                                  instrument=instName,
+                                  physical_filter=physical_filter)
+            tputs = passband.sampleAt(position=geom.Point2D(0.0, 0.0), wavelengths=lambdas)
+            self.assertEqual(np.min(tputs), 0.0)
+            self.assertGreater(np.max(tputs), 0.0)
+
     def _testFgcmBuildStarsTable(self, instName, testName, queryString, visits, nStar, nObs):
         """Test running of FgcmBuildStarsTableTask