Merge pull request #15 from lsst/tickets/DM-6622

Review complete. Tickets/dm 6622

include/lsst/jointcal/Jointcal.h

@@ -20,8 +20,9 @@ namespace jointcal {
     struct JointcalControl {
         LSST_CONTROL_FIELD(sourceFluxField, std::string, "name of flux field in source catalog");
 
-        JointcalControl() :
-            sourceFluxField("None")
+        JointcalControl(std::string const sourceFluxType) :
+            // Set sourceFluxType to the value used in the source selector.
+            sourceFluxField("slot_"+sourceFluxType+"Flux")
         {
             validate();
         }

python/lsst/jointcal/jointcal.py

@@ -89,12 +89,6 @@ class JointcalConfig(pexConfig.Config):
         dtype = int,
         default = 3,
     )
-    sourceFluxType = pexConfig.Field(
-        doc = "Type of source flux (e.g. Ap, Psf, Calib): passed to sourceSelector "
-              "and used in ccdImage",
-        dtype = str,
-        default = "Calib"
-    )
     sourceSelector = sourceSelectorRegistry.makeField(
         doc = "How to select sources for cross-matching",
         default = "astrometry"
@@ -103,9 +97,10 @@ class JointcalConfig(pexConfig.Config):
     def setDefaults(self):
         sourceSelector = self.sourceSelector["astrometry"]
         sourceSelector.setDefaults()
-        sourceSelector.sourceFluxType = self.sourceFluxType
         # don't want to lose existing flags, just add to them.
         sourceSelector.badFlags.extend(["slot_Shape_flag"])
+        # This should be used to set the FluxField value in jointcal::JointcalControl
+        sourceSelector.sourceFluxType = 'Calib'
 
 
 class JointcalTask(pipeBase.CmdLineTask):
@@ -213,9 +208,7 @@ def run(self, dataRefs, profile_jointcal=False):
         if len(dataRefs) == 0:
             raise ValueError('Need a list of data references!')
 
-        jointcalControl = jointcalLib.JointcalControl()
-        jointcalControl.sourceFluxField = 'slot_'+self.config.sourceFluxType+'Flux'
-
+        jointcalControl = jointcalLib.JointcalControl(self.sourceSelector.config.sourceFluxType)
         associations = jointcalLib.Associations()
 
         load_cat_prof_file = 'jointcal_load_catalog.prof' if profile_jointcal else ''
@@ -253,7 +246,6 @@ def run(self, dataRefs, profile_jointcal=False):
         # Determine default filter associated to the catalog
         filt, mfilt = andConfig.magColumnMap.items()[0]
         print("Using", filt, "band for reference flux")
-
         refCat = loader.loadSkyCircle(center, afwGeom.Angle(radius, afwGeom.radians), filt).refCat
 
         # associations.CollectRefStars(False) # To use USNO-A catalog

src/CcdImage.cc

@@ -18,38 +18,17 @@ namespace afwImg = lsst::afw::image;
 namespace lsst {
 namespace jointcal {
 
-  /* Interesting fields of the stack catalogs :
- 'base_SdssCentroid_x'
- 'base_SdssCentroid_y'
- 'base_SdssCentroid_xSigma'
- 'base_SdssCentroid_ySigma'
-
-  We miss the xy uncertainty term.
-  We can cook it up from the sdss shape:
-  'base_SdssShape_xx'
-  'base_SdssShape_yy'
-  'base_SdssShape_xy'
-
-  for fluxes, we might use :
-  'base_CircularApertureFlux_2_flux'
-  'base_CircularApertureFlux_2_fluxSigma'
-
-   where the '2' should be read from the environment.
-  */
-
 static double sq(const double &x) { return x*x;}
 
 void CcdImage::LoadCatalog(const lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceRecord> &Cat,const std::string &fluxField)
 {
-  // TODO: these should use slots, e.g. "slot_Centroid_x", or better:
-  // self.config.centroidName + "_x" in python.
-  auto xKey = Cat.getSchema().find<double>("base_SdssCentroid_x").key;
-  auto yKey = Cat.getSchema().find<double>("base_SdssCentroid_y").key;
-  auto xsKey = Cat.getSchema().find<float>("base_SdssCentroid_xSigma").key;
-  auto ysKey = Cat.getSchema().find<float>("base_SdssCentroid_ySigma").key;
-  auto mxxKey = Cat.getSchema().find<double>("base_SdssShape_xx").key;
-  auto myyKey = Cat.getSchema().find<double>("base_SdssShape_yy").key;
-  auto mxyKey = Cat.getSchema().find<double>("base_SdssShape_xy").key;
+  auto xKey = Cat.getSchema().find<double>("slot_Centroid_x").key;
+  auto yKey = Cat.getSchema().find<double>("slot_Centroid_y").key;
+  auto xsKey = Cat.getSchema().find<float>("slot_Centroid_xSigma").key;
+  auto ysKey = Cat.getSchema().find<float>("slot_Centroid_ySigma").key;
+  auto mxxKey = Cat.getSchema().find<double>("slot_Shape_xx").key;
+  auto myyKey = Cat.getSchema().find<double>("slot_Shape_yy").key;
+  auto mxyKey = Cat.getSchema().find<double>("slot_Shape_xy").key;
   auto fluxKey = Cat.getSchema().find<double>(fluxField + "_flux").key;
   auto efluxKey = Cat.getSchema().find<double>(fluxField  + "_fluxSigma").key;
 
@@ -188,7 +167,7 @@ CcdImage::CcdImage(lsst::afw::table::SortedCatalogT<lsst::afw::table::SourceReco
         ra = lsst::afw::geom::degToRad(meta->get<double>("RA_DEG"));
         dec = lsst::afw::geom::degToRad(meta->get<double>("DEC_DEG"));
     }
-    else if (camera == "hsc") {
+    else if (camera == "HSC") {
         airMass = meta->get<double>("AIRMASS");
         jd = meta->get<double>("MJD");  // Julian date
         expTime = meta->get<double>("EXPTIME");

tests/testJointcal_hsc.py

@@ -0,0 +1,74 @@
+# See COPYRIGHT file at the top of the source tree.
+
+from __future__ import division, absolute_import, print_function
+
+import unittest
+import os
+
+from astropy import units as u
+
+import lsst.afw.geom
+import lsst.afw.coord as afwCoord
+import lsst.utils
+import lsst.pex.exceptions
+
+import jointcalTestBase
+
+try:
+    data_dir = lsst.utils.getPackageDir('validation_data_hsc')
+    os.environ['ASTROMETRY_NET_DATA_DIR'] = os.path.join(data_dir, 'sdss-dr9-fink-v5b')
+except lsst.pex.exceptions.NotFoundError:
+    data_dir = None
+
+# We don't want the absolute astrometry to become significantly worse
+# than the single-epoch astrometry (about 0.040").
+# This value was empirically determined from the first run of jointcal on
+# this data, and will likely vary from survey to survey.
+absolute_error = 52e-3*u.arcsecond
+# Set to True for a comparison plot and some diagnostic numbers.
+do_plot = False
+
+
+# for MemoryTestCase
+def setup_module(module):
+    lsst.utils.tests.init()
+
+
+class JointcalTestHSC(jointcalTestBase.JointcalTestBase, lsst.utils.tests.TestCase):
+    def setUp(self):
+        jointcalTestBase.JointcalTestBase.setUp(self)
+        self.do_plot = do_plot
+        self.match_radius = 0.1*lsst.afw.geom.arcseconds
+
+        # position of this validation_data_hsc catalog
+        center = afwCoord.IcrsCoord(320.367492*lsst.afw.geom.degrees, 0.3131554*lsst.afw.geom.degrees)
+        radius = 5*lsst.afw.geom.degrees
+        self._prep_reference_loader(center, radius)
+
+        self.input_dir = os.path.join(data_dir, 'DATA')
+        self.all_visits = [903334, 903336, 903338, 903342, 903344, 903346]
+
+    @unittest.skipIf(data_dir is None, "validation_data_hsc not setup")
+    def test_jointcalTask_2_visits(self):
+        # NOTE: The relative RMS limit was empirically determined from the
+        # first run of jointcal on this data. We should always do better than
+        # this in the future!
+        relative_error = 17e-3*u.arcsecond
+        self._testJointCalTask(2, relative_error, absolute_error)
+
+    @unittest.skipIf(data_dir is None, "validation_data_hsc not setup")
+    def test_jointcalTask_6_visits(self):
+        # NOTE: The relative RMS limit was empirically determined from the
+        # first run of jointcal on this data. We should always do better than
+        # this in the future!
+        relative_error = 10e-3*u.arcsecond
+        self._testJointCalTask(6, relative_error, absolute_error)
+
+
+# TODO: the memory test cases currently fail in jointcal. Filed as DM-6626.
+# class MyMemoryTestCase(lsst.utils.tests.MemoryTestCase):
+#     pass
+
+if __name__ == "__main__":
+    lsst.utils.tests.init()
+    unittest.main()

tests/testJointcal_lsstSim.py

@@ -106,7 +106,7 @@ def setUp(self):
         # sources used for associations, and astrometrySourceSelector does not
         # exactly match the original bundled StarSelector.
         # TODO: Once we make jointcal more robust, we should be able to drop this.
-        self.jointcalTask.config.sourceSelector["astrometry"].minSnr = 13
+        self.jointcalTask.sourceSelector.config.minSnr = 13
 
     @unittest.skipIf(data_dir is None, "validation_data_jointcal not setup")
     @unittest.skip('jointcal currently fails if only given one catalog!')

ups/jointcal.table

@@ -13,6 +13,7 @@ setupRequired(jointcal_cholmod)
 
 # for the testing files
 setupOptional(validation_data_jointcal)
+setupOptional(validation_data_hsc)
 
 envPrepend(LD_LIBRARY_PATH, ${PRODUCT_DIR}/lib)
 envPrepend(DYLD_LIBRARY_PATH, ${PRODUCT_DIR}/lib)