Set magnitude zero point if none exists

This code will now check that calib objects have a valid fluxMag0 (zero-point) set. If not a default calib object is used with a reference zero-point of 27.
lsst · May 9, 2016 · 8abe2f1 · 8abe2f1
1 parent ac94e06
commit 8abe2f1
Show file tree

Hide file tree

Showing 3 changed files with 45 additions and 7 deletions.
diff --git a/include/lsst/meas/modelfit/UnitSystem.h b/include/lsst/meas/modelfit/UnitSystem.h
@@ -24,6 +24,8 @@
 #ifndef LSST_MEAS_MODELFIT_UnitSystem_h_INCLUDED
 #define LSST_MEAS_MODELFIT_UnitSystem_h_INCLUDED
 
+#include <math.h>
+
 #include "lsst/afw/image/Exposure.h"
 #include "lsst/afw/image/Wcs.h"
 #include "lsst/afw/image/Calib.h"
@@ -48,11 +50,13 @@ struct UnitSystem {
      *  set such that unit flux is the given magnitude.  See @ref modelfitUnits for an explanation
      *  of why we frequently use this system.
      */
+    UnitSystem(afw::coord::Coord const & position, std::shared_ptr<const lsst::afw::image::Calib> calibIn,
+               double flux);
     UnitSystem(afw::coord::Coord const & position, Scalar mag);
 
     /// Construct a UnitSystem from a give Wcs and Calib
     UnitSystem(PTR(afw::image::Wcs const) wcs_, PTR(afw::image::Calib const) calib_) :
-        wcs(wcs_), calib(calib_)
+        wcs(wcs_), calib(validateCalib(calib_))
     {}
 
 // work around a bug in SWIG 3.0.2: mis-handling templated constructors
@@ -61,17 +65,21 @@ struct UnitSystem {
     /// Construct a UnitSystem by extracting the Wcs and Calib from an Exposure (implicit)
     template <typename T>
     UnitSystem(afw::image::Exposure<T> const & exposure) :
-        wcs(exposure.getWcs()), calib(exposure.getCalib())
+        wcs(exposure.getWcs()), calib(validateCalib(exposure.getCalib()))
     {}
 #else
     UnitSystem(afw::image::Exposure<float> const & exposure) :
-        wcs(exposure.getWcs()), calib(exposure.getCalib())
+        wcs(exposure.getWcs()), calib(validateCalib(exposure.getCalib()))
     {}
     UnitSystem(afw::image::Exposure<double> const & exposure) :
-        wcs(exposure.getWcs()), calib(exposure.getCalib())
+        wcs(exposure.getWcs()), calib(validateCalib(exposure.getCalib()))
     {}
-
 #endif
+
+private:
+    std::shared_ptr<const lsst::afw::image::Calib> validateCalib(
+                                                    std::shared_ptr<const lsst::afw::image::Calib> calib_);
+    static std::shared_ptr<const lsst::afw::image::Calib> getDefaultCalib();
 };
 
 /**

diff --git a/src/CModel.cc b/src/CModel.cc
@@ -555,7 +555,7 @@ struct CModelStageData {
         Model const & model
     ) :
         measSysCenter(center), position(exposure.getWcs()->pixelToSky(center)),
-        measSys(exposure), fitSys(*position, exposure.getCalib()->getMagnitude(approxFlux)),
+        measSys(exposure), fitSys(*position, exposure.getCalib(), approxFlux),
         fitSysToMeasSys(*position, fitSys, measSys),
         parameters(ndarray::allocate(model.getNonlinearDim() + model.getAmplitudeDim())),
         nonlinear(parameters[ndarray::view(0, model.getNonlinearDim())]),

diff --git a/src/UnitSystem.cc b/src/UnitSystem.cc
@@ -28,14 +28,44 @@
 
 namespace lsst { namespace meas { namespace modelfit {
 
-UnitSystem::UnitSystem(afw::coord::Coord const & position, Scalar mag) {
+UnitSystem::UnitSystem(afw::coord::Coord const & position,
+                       std::shared_ptr<const lsst::afw::image::Calib> calibIn,
+                       double flux) {
     double cdelt = (1.0*lsst::afw::geom::arcseconds).asDegrees();
     wcs = afw::image::makeWcs(position, afw::geom::Point2D(0.0, 0.0), cdelt, 0.0, 0.0, cdelt);
+    calibIn = validateCalib(calibIn);
+    Scalar mag = calibIn->getMagnitude(flux);
     PTR(afw::image::Calib) calib_ = std::make_shared<afw::image::Calib>();
     calib_->setFluxMag0(std::pow(10.0, mag/2.5));
     calib = calib_;
 }
 
+UnitSystem::UnitSystem(afw::coord::Coord const & position, Scalar mag) {
+        double cdelt = (1.0*lsst::afw::geom::arcseconds).asDegrees();
+        wcs = afw::image::makeWcs(position, afw::geom::Point2D(0.0, 0.0), cdelt, 0.0, 0.0, cdelt);
+        PTR(afw::image::Calib) calib_ = std::make_shared<afw::image::Calib>();
+        calib_->setFluxMag0(std::pow(10.0, mag/2.5));
+        calib = calib_;
+    }
+
+std::shared_ptr<const lsst::afw::image::Calib> UnitSystem::validateCalib(
+        std::shared_ptr<const lsst::afw::image::Calib> calib_) {
+        if (calib_->getFluxMag0().first == 0.0){
+            return getDefaultCalib();
+        } else{
+            return calib_;
+        }
+    }
+
+std::shared_ptr<const lsst::afw::image::Calib> UnitSystem::getDefaultCalib() {
+        // Create a default calib object with a zero-point set to magnitude 27
+        static std::shared_ptr<const lsst::afw::image::Calib> tmp =
+                                                            std::make_shared<const lsst::afw::image::Calib>(
+                                                            std::pow(10.0, 27.0/2.5)
+                                                            );
+        return tmp;
+    }
+
 LocalUnitTransform::LocalUnitTransform(
     afw::coord::Coord const & position,
     UnitSystem const & source,