Merge pull request #48 from lsst/tickets/DM-16650

DM-16650: Update to new nJy PhotoCalib definition

python/lsst/meas/mosaic/mosaicTask.py

@@ -701,12 +701,14 @@ def writeFcr(self, dataRefList):
                 print("failed to read Wcs for PhotoCalib: %s" % (e))
                 continue
             instFluxMag0, instFluxMag0Err = calib.getFluxMag0()
+            # TODO: need to scale these until DM-10153 is completed and PhotoCalib has replaced Calib
+            referenceFlux = 1e23 * 10**(48.6 / -2.5) * 1e9
             bf = measMosaic.FluxFitBoundedField(bbox, newP, wcs,
                                                 zeroPoint=instFluxMag0,
                                                 nQuarter=nQuarter)
             photoCalib = afwImage.PhotoCalib(
-                1.0/instFluxMag0,
-                instFluxMag0Err/instFluxMag0**2,
+                referenceFlux / instFluxMag0,
+                referenceFlux * instFluxMag0Err/instFluxMag0**2,
                 bf,
                 isConstant=False
             )

src/FluxFitBoundedField.cc

@@ -26,6 +26,7 @@
 #include "lsst/afw/table/io/InputArchive.h"
 #include "lsst/afw/table/io/CatalogVector.h"
 #include "lsst/afw/table/io/Persistable.cc"
+#include "lsst/utils/Magnitude.h"
 
 namespace lsst {
 namespace afw {
@@ -89,7 +90,7 @@ FluxFitBoundedField::FluxFitBoundedField(
 }
 
 double FluxFitBoundedField::evaluate(afw::geom::Point2D const & position) const {
-    double r = 1.0/_zeroPoint;
+    double r = utils::referenceFlux/_zeroPoint;
     if (_ffp) {
         auto xy = _transform(position);
         if (xy.getX() < -1E-8 || xy.getY() < -1E-8) {

tests/test_fluxFitBoundedField.py

@@ -36,6 +36,10 @@
 
 DATA_DIR = os.path.join(os.path.split(__file__)[0], "data")
 
+# Note: mosaic's internal ffp object is magnitude-based.
+# Instead, for these tests, we just divide this out when comparing.
+referenceFlux = 1e23 * 10**(48.6 / -2.5) * 1e9
+
 
 def displayImageDifferences(image1, image2, rtol=1E-8, atol=1E-8, pause=False):
     import lsst.afw.display
@@ -119,7 +123,7 @@ class FluxFitBoundedFieldTestCase(lsst.utils.tests.TestCase):
     def setUp(self):
         np.random.seed(100)
         self.tract = 8766
-        self.visit = 7358
+        self.visit = 11506
         self.ccds = {0: 49, 2: 50, 3: 101, 1: 102}
         # Box is the same for all CCDs, since it's defined in CCD coordinates,
         # which are rotated w.r.t. focal plane coordinates.
@@ -142,21 +146,18 @@ def setUp(self):
             self.ffp[ccd] = lsst.meas.mosaic.FluxFitParams(fcrMetadata)
             wcsFilename = os.path.join(
                 DATA_DIR,
-                "%d/wcs-%07d-%03d.fits" % (self.tract, self.visit, ccd)
+                "%d/jointcal_wcs-%07d-%03d.fits" % (self.tract, self.visit, ccd)
             )
-            # WCS is saved as an empty exposure with WCS
-            # (rather than as a WCS fits file)
-            wcsExposure = lsst.afw.image.ExposureF(wcsFilename)
-            self.wcs[ccd] = wcsExposure.getWcs()
+            self.wcs[ccd] = lsst.afw.geom.SkyWcs.readFits(wcsFilename)
             photoCalibFilename = os.path.join(
                 DATA_DIR,
-                "%d/photoCalib-%07d-%03d.fits" % (self.tract, self.visit, ccd)
+                "%d/jointcal_photoCalib-%07d-%03d.fits" % (self.tract, self.visit, ccd)
             )
             self.photoCalib[ccd] = lsst.afw.image.PhotoCalib.readFits(photoCalibFilename)
             camera[ccd] = MockDetector(MockOrientation(nQuarter))
             self.dataRefs[ccd] = MockDataRef(visit=self.visit, tract=self.tract, ccd=ccd)
             self.dataRefs[ccd].put(fcrMetadata, "fcr_md", )
-            self.dataRefs[ccd].put(wcsExposure, "wcs")
+            self.dataRefs[ccd].put(self.wcs[ccd], "jointcal_wcs")
             self.dataRefs[ccd].put(calexpMetadata, "calexp_md")
             self.dataRefs[ccd].put(camera, "camera")
             self.dataRefs[ccd].put(self.bbox, "calexp_bbox")
@@ -189,6 +190,7 @@ def checkFillImage(self, nQuarter, ffp=True, wcs=True, display=False):
         bf, ffp, wcs = self.makeBoundedField(nQuarter, ffp, wcs)
         image1 = lsst.afw.image.ImageF(self.bbox)
         bf.fillImage(image1, xStep=100, yStep=16)
+        image1 /= referenceFlux  # turn nJy into Maggy: ffp outputs Maggy
         if nQuarter%2:
             width, height = self.bbox.getHeight(), self.bbox.getWidth()
         else:
@@ -265,9 +267,9 @@ def checkPhotoCalibCatalog(self, nQuarter):
         mag2 = catalog2["example_mag"]
         # Check that the non-spatially varying part of the correction is the same.
         fluxMag0 = results2.ffp.calib.getFluxMag0()
-        self.assertFloatsAlmostEqual(photoCalib.getInstFluxMag0(), fluxMag0[0],
+        self.assertFloatsAlmostEqual(photoCalib.getInstFluxAtZeroMagnitude(), fluxMag0[0],
                                      rtol=1E-14)
-        self.assertFloatsAlmostEqual(photoCalib.getCalibrationErr(), fluxMag0[1]/fluxMag0[0]**2,
+        self.assertFloatsAlmostEqual(photoCalib.getCalibrationErr(), referenceFlux*fluxMag0[1]/fluxMag0[0]**2,
                                      rtol=1E-14)
         # Compute partially-calibrated magnitudes that don't account for the spatially-varying part.
         mag0, magErr0 = results2.ffp.calib.getMagnitude(catalog.get("example_instFlux"),
@@ -279,7 +281,7 @@ def checkPhotoCalibCatalog(self, nQuarter):
         # round-off error).
         magDiff2 = mag2 - mag0
         magDiff1 = mag1 - mag0
-        self.assertFloatsAlmostEqual(magDiff1, magDiff2, rtol=rtol*2E3)
+        self.assertFloatsAlmostEqual(magDiff1, magDiff2, rtol=rtol*3E3)
 
     def checkPhotoCalibExposure(self, nQuarter):
         ccd = self.ccds[nQuarter]