LSSTApplications  10.0-2-g4f67435,11.0.rc2+1,11.0.rc2+12,11.0.rc2+3,11.0.rc2+4,11.0.rc2+5,11.0.rc2+6,11.0.rc2+7,11.0.rc2+8
LSSTDataManagementBasePackage
Calib.cc
Go to the documentation of this file.
1 // -*- lsst-c++ -*-
2 
3 /*
4  * LSST Data Management System
5  * Copyright 2008, 2009, 2010 LSST Corporation.
6  *
7  * This product includes software developed by the
8  * LSST Project (http://www.lsst.org/).
9  *
10  * This program is free software: you can redistribute it and/or modify
11  * it under the terms of the GNU General Public License as published by
12  * the Free Software Foundation, either version 3 of the License, or
13  * (at your option) any later version.
14  *
15  * This program is distributed in the hope that it will be useful,
16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18  * GNU General Public License for more details.
19  *
20  * You should have received a copy of the LSST License Statement and
21  * the GNU General Public License along with this program. If not,
22  * see <http://www.lsstcorp.org/LegalNotices/>.
23  */
24 
30 #include <cmath>
31 
32 #include "boost/format.hpp"
33 #include "boost/lexical_cast.hpp"
34 #include "boost/algorithm/string/trim.hpp"
35 
36 #include "ndarray.h"
37 #include "lsst/pex/exceptions.h"
39 #include "lsst/afw/image/Calib.h"
44 
45 namespace lsst { namespace afw { namespace image {
49 Calib::Calib() : _midTime(), _exptime(0.0), _fluxMag0(0.0), _fluxMag0Sigma(0.0) {}
55 Calib::Calib(std::vector<CONST_PTR(Calib)> const& calibs
56  ) :
57  _midTime(), _exptime(0.0), _fluxMag0(0.0), _fluxMag0Sigma(0.0)
58 {
59  if (calibs.empty()) {
60  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
61  "You must provide at least one input Calib");
62  }
63 
64  double const fluxMag00 = calibs[0]->_fluxMag0;
65  double const fluxMag0Sigma0 = calibs[0]->_fluxMag0Sigma;
66 
67  double midTimeSum = 0.0; // sum(time*expTime)
68  for (std::vector<CONST_PTR(Calib)>::const_iterator ptr = calibs.begin(); ptr != calibs.end(); ++ptr) {
69  Calib const& calib = **ptr;
70 
71  if (::fabs(fluxMag00 - calib._fluxMag0) > std::numeric_limits<double>::epsilon() ||
72  ::fabs(fluxMag0Sigma0 - calib._fluxMag0Sigma) > std::numeric_limits<double>::epsilon()) {
73  throw LSST_EXCEPT(lsst::pex::exceptions::InvalidParameterError,
74  (boost::format("You may only combine calibs with the same fluxMag0: "
75  "%g +- %g v %g +- %g")
76  % calib.getFluxMag0().first % calib.getFluxMag0().second
77  % calibs[0]->getFluxMag0().first % calibs[0]->getFluxMag0().second
78  ).str());
79  }
80 
81  double const exptime = calib._exptime;
82 
83  midTimeSum += calib._midTime.get()*exptime;
84  _exptime += exptime;
85  }
86 
87  daf::base::DateTime tmp(midTimeSum/_exptime); // there's no way to set the double value directly
88  using std::swap;
89  swap(_midTime, tmp);
90 }
91 
97  double exptime = 0.0;
98  double fluxMag0 = 0.0, fluxMag0Sigma = 0.0;
99 
100  std::string key = "TIME-MID";
101  if (metadata->exists(key)) {
102  midTime = lsst::daf::base::DateTime(boost::algorithm::trim_right_copy(metadata->getAsString(key)));
103  }
104 
105  key = "EXPTIME";
106  if (metadata->exists(key)) {
107  try {
108  exptime = metadata->getAsDouble(key);
109  } catch (lsst::pex::exceptions::TypeError & err) {
110  std::string exptimeStr = metadata->getAsString(key);
111  exptime = boost::lexical_cast<double>(exptimeStr);
112  }
113  }
114 
115  key = "FLUXMAG0";
116  if (metadata->exists(key)) {
117  fluxMag0 = metadata->getAsDouble(key);
118 
119  key = "FLUXMAG0ERR";
120  if (metadata->exists(key)) {
121  fluxMag0Sigma = metadata->getAsDouble(key);
122  }
123  }
124 
125  _midTime = midTime;
126  _exptime = exptime;
129 }
133 bool Calib::_throwOnNegativeFlux = true;
137 void
138 Calib::setThrowOnNegativeFlux(bool raiseException
139  )
140 {
141  _throwOnNegativeFlux = raiseException;
142 }
143 
147 bool
149 {
150  return _throwOnNegativeFlux;
151 }
152 
153 namespace detail {
160  )
161 {
162  int nstripped = 0;
163 
164  std::string key = "TIME-MID";
165  if (metadata->exists(key)) {
166  metadata->remove(key);
167  nstripped++;
168  }
169 
170  key = "EXPTIME";
171  if (metadata->exists(key)) {
172  metadata->remove(key);
173  nstripped++;
174  }
175 
176  key = "FLUXMAG0";
177  if (metadata->exists(key)) {
178  metadata->remove(key);
179  nstripped++;
180  }
181 
182  key = "FLUXMAG0ERR";
183  if (metadata->exists(key)) {
184  metadata->remove(key);
185  nstripped++;
186  }
187 
188  return nstripped;
189 }
190 }
191 
197 bool Calib::operator==(Calib const& rhs) const {
198  return
199  ::fabs(_midTime.get() - rhs._midTime.get()) < std::numeric_limits<double>::epsilon() &&
200  _exptime == rhs._exptime &&
201  _fluxMag0 == rhs._fluxMag0 &&
203 }
204 
211  )
212 {
213  _midTime = midTime;
214 }
215 
220 {
221  return _midTime;
222 }
223 
231  boost::shared_ptr<const lsst::afw::cameraGeom::Detector>,
233  ) const
234 {
235  return _midTime;
236 }
237 
241 void Calib::setExptime(double exptime
242  ) {
243  _exptime = exptime;
244 }
245 
249 double Calib::getExptime() const
250 {
251  return _exptime;
252 }
253 
258  double fluxMag0Sigma
259  )
260 {
263 }
264 void Calib::setFluxMag0(std::pair<double, double> fluxMag0AndSigma
265  )
266 {
267  _fluxMag0 = fluxMag0AndSigma.first;
268  _fluxMag0Sigma = fluxMag0AndSigma.second;
269 }
270 
274 std::pair<double, double> Calib::getFluxMag0() const
275 {
276  return std::make_pair(_fluxMag0, _fluxMag0Sigma);
277 }
278 
279 namespace {
280 inline void checkNegativeFlux0(double fluxMag0) {
281  if (fluxMag0 <= 0) {
282  throw LSST_EXCEPT(lsst::pex::exceptions::DomainError,
283  (boost::format("Flux of 0-mag object must be >= 0: saw %g") % fluxMag0).str());
284  }
285 }
286 inline bool isNegativeFlux(double flux, bool doThrow)
287 {
288  if (flux <= 0) {
289  if (doThrow) {
290  throw LSST_EXCEPT(lsst::pex::exceptions::DomainError,
291  (boost::format("Flux must be >= 0: saw %g") % flux).str());
292  }
293  return true;
294  }
295  return false;
296 }
297 inline double convertToFlux(double fluxMag0, double mag) {
298  return fluxMag0 * ::pow(10.0, -0.4 * mag);
299 }
300 inline double convertToFluxErr(double fluxMag0InvSNR, double flux, double magErr) {
301  // Want to:
302  // return flux * hypot(_fluxMag0Sigma/_fluxMag0, 0.4*std::log(10)*magSigma/mag);
303  // But hypot is not standard C++ so use <http://en.wikipedia.org/wiki/Hypot#Implementation>
304  double a = fluxMag0InvSNR;
305  double b = 0.4 * std::log(10.0) * magErr;
306  if (std::abs(a) < std::abs(b)) {
307  std::swap(a, b);
308  }
309  return flux * std::abs(a) * std::sqrt(1 + std::pow(b / a, 2));
310 }
311 inline double convertToMag(double fluxMag0, double flux) {
312  return -2.5*::log10(flux/fluxMag0);
313 }
314 inline void convertToMagWithErr(double *mag, double *magErr, double fluxMag0, double fluxMag0InvSNR,
315  double flux, double fluxErr)
316 {
317  double const rat = flux/fluxMag0;
318  double const ratErr = ::sqrt((::pow(fluxErr, 2) + ::pow(flux*fluxMag0InvSNR, 2))/::pow(fluxMag0, 2));
319 
320  *mag = -2.5*::log10(rat);
321  *magErr = 2.5/::log(10.0)*ratErr/rat;
322 }
323 
324 } // anonymous namespace
325 
329 double Calib::getFlux(double const mag
330  ) const {
331  checkNegativeFlux0(_fluxMag0);
332  return convertToFlux(_fluxMag0, mag);
333 }
335  checkNegativeFlux0(_fluxMag0);
338  ndarray::Array<double,1>::Iterator outIter = flux.begin();
339  for (; inIter != mag.end(); ++inIter, ++outIter) {
340  *outIter = convertToFlux(_fluxMag0, *inIter);
341  }
342  return flux;
343 }
344 
350 std::pair<double, double> Calib::getFlux(
351  double const mag,
352  double const magSigma
353  ) const
354 {
355  checkNegativeFlux0(_fluxMag0);
356  double const flux = convertToFlux(_fluxMag0, mag);
357  double const fluxErr = convertToFluxErr(_fluxMag0Sigma/_fluxMag0, flux, magSigma);
358  return std::make_pair(flux, fluxErr);
359 }
360 
361 std::pair<ndarray::Array<double,1>, ndarray::Array<double,1> > Calib::getFlux(
362  ndarray::Array<double const,1> const & mag,
363  ndarray::Array<double const,1> const & magErr
364 ) const {
365  checkNegativeFlux0(_fluxMag0);
366  if (mag.size() != magErr.size()) {
367  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
368  (boost::format("Size of mag (%d) and magErr (%d) don't match") %
369  mag.size() % magErr.size()).str());
370  }
371 
375  ndarray::Array<double const,1>::Iterator magErrIter = magErr.begin();
376  ndarray::Array<double,1>::Iterator fluxIter = flux.begin();
377  ndarray::Array<double,1>::Iterator fluxErrIter = fluxErr.begin();
378 
379  double fluxMag0InvSNR = _fluxMag0Sigma/_fluxMag0;
380  for (; magIter != mag.end(); ++magIter, ++magErrIter, ++fluxIter, ++fluxErrIter) {
381  *fluxIter = convertToFlux(_fluxMag0, *magIter);
382  *fluxErrIter = convertToFluxErr(fluxMag0InvSNR, *fluxIter, *magErrIter);
383  }
384 
385  return std::make_pair(flux, fluxErr);
386 }
387 
391 double Calib::getMagnitude(double const flux
392  ) const
393 {
394  checkNegativeFlux0(_fluxMag0);
395  if (isNegativeFlux(flux, Calib::getThrowOnNegativeFlux())) {
396  return std::numeric_limits<double>::quiet_NaN();
397  }
398  return convertToMag(_fluxMag0, flux);
399 }
400 
404 std::pair<double, double> Calib::getMagnitude(double const flux,
405  double const fluxErr
406  ) const
407 {
408  checkNegativeFlux0(_fluxMag0);
409  if (isNegativeFlux(flux, Calib::getThrowOnNegativeFlux())) {
410  double const NaN = std::numeric_limits<double>::quiet_NaN();
411  return std::make_pair(NaN, NaN);
412  }
413 
414  double mag, magErr;
415  convertToMagWithErr(&mag, &magErr, _fluxMag0, _fluxMag0Sigma/_fluxMag0, flux, fluxErr);
416  return std::make_pair(mag, magErr);
417 }
418 
420  checkNegativeFlux0(_fluxMag0);
424  int nonPositive = 0;
425  for (; fluxIter != flux.end(); ++fluxIter, ++magIter) {
426  if (isNegativeFlux(*fluxIter, false)) {
427  ++nonPositive;
428  *magIter = std::numeric_limits<double>::quiet_NaN();
429  continue;
430  }
431  *magIter = convertToMag(_fluxMag0, *fluxIter);
432  }
433  if (nonPositive && Calib::getThrowOnNegativeFlux()) {
434  throw LSST_EXCEPT(lsst::pex::exceptions::DomainError,
435  (boost::format("Flux must be >= 0: %d non-positive seen") % nonPositive).str());
436  }
437  return mag;
438 }
439 
440 std::pair<ndarray::Array<double,1>, ndarray::Array<double,1> > Calib::getMagnitude(
441  ndarray::Array<double const,1> const & flux,
442  ndarray::Array<double const,1> const & fluxErr
443 ) const {
444  checkNegativeFlux0(_fluxMag0);
445  if (flux.size() != fluxErr.size()) {
446  throw LSST_EXCEPT(lsst::pex::exceptions::LengthError,
447  (boost::format("Size of flux (%d) and fluxErr (%d) don't match") %
448  flux.size() % fluxErr.size()).str());
449  }
450 
454  ndarray::Array<double const,1>::Iterator fluxErrIter = fluxErr.begin();
456  ndarray::Array<double,1>::Iterator magErrIter = magErr.begin();
457  int nonPositive = 0;
458  double fluxMag0InvSNR = _fluxMag0Sigma/_fluxMag0;
459  for (; fluxIter != flux.end(); ++fluxIter, ++fluxErrIter, ++magIter, ++magErrIter) {
460  if (isNegativeFlux(*fluxIter, false)) {
461  ++nonPositive;
462  double const NaN = std::numeric_limits<double>::quiet_NaN();
463  *magIter = NaN;
464  *magErrIter = NaN;
465  continue;
466  }
467  double f, df;
468  convertToMagWithErr(&f, &df, _fluxMag0, fluxMag0InvSNR, *fluxIter, *fluxErrIter);
469  *magIter = f;
470  *magErrIter = df;
471  }
472  if (nonPositive && Calib::getThrowOnNegativeFlux()) {
473  throw LSST_EXCEPT(lsst::pex::exceptions::DomainError,
474  (boost::format("Flux must be >= 0: %d non-positive seen") % nonPositive).str());
475  }
476  return std::make_pair(mag, magErr);
477 }
478 
479 namespace {
480 
481 class CalibSchema : private boost::noncopyable {
482 public:
483  table::Schema schema;
484  table::Key<boost::int64_t> midTime;
485  table::Key<double> expTime;
486  table::Key<double> fluxMag0;
487  table::Key<double> fluxMag0Sigma;
488 
489  static CalibSchema const & get() {
490  static CalibSchema instance;
491  return instance;
492  }
493 
494 private:
495  CalibSchema() :
496  schema(),
497  midTime(schema.addField<boost::int64_t>(
498  "midtime", "middle of the time of the exposure relative to Unix epoch", "nanoseconds"
499  )),
500  expTime(schema.addField<double>("exptime", "exposure time", "seconds")),
501  fluxMag0(schema.addField<double>("fluxmag0", "flux of a zero-magnitude object", "dn")),
502  fluxMag0Sigma(schema.addField<double>("fluxmag0.err", "1-sigma error on fluxmag0", "dn"))
503  {
504  schema.getCitizen().markPersistent();
505  }
506 };
507 
508 class CalibFactory : public table::io::PersistableFactory {
509 public:
510 
511  virtual PTR(table::io::Persistable)
512  read(InputArchive const & archive, CatalogVector const & catalogs) const {
513  CalibSchema const & keys = CalibSchema::get();
514  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
515  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
516  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
517  table::BaseRecord const & record = catalogs.front().front();
518  PTR(Calib) result(new Calib());
519  result->setMidTime(daf::base::DateTime(static_cast<long long>(record.get(keys.midTime))));
520  result->setExptime(record.get(keys.expTime));
521  result->setFluxMag0(record.get(keys.fluxMag0), record.get(keys.fluxMag0Sigma));
522  return result;
523  }
524 
525  explicit CalibFactory(std::string const & name) : table::io::PersistableFactory(name) {}
526 
527 };
528 
529 std::string getCalibPersistenceName() { return "Calib"; }
530 
531 CalibFactory registration(getCalibPersistenceName());
532 
533 } // anonymous
534 
535 std::string Calib::getPersistenceName() const { return getCalibPersistenceName(); }
536 
537 void Calib::write(OutputArchiveHandle & handle) const {
538  CalibSchema const & keys = CalibSchema::get();
539  table::BaseCatalog cat = handle.makeCatalog(keys.schema);
540  PTR(table::BaseRecord) record = cat.addNew();
541  record->set(keys.midTime, getMidTime().nsecs());
542  record->set(keys.expTime, getExptime());
543  std::pair<double,double> fluxMag0 = getFluxMag0();
544  record->set(keys.fluxMag0, fluxMag0.first);
545  record->set(keys.fluxMag0Sigma, fluxMag0.second);
546  handle.saveCatalog(cat);
547 }
548 
549 }}} // lsst::afw::image
double getMagnitude(double const flux) const
Definition: Calib.cc:391
Class for handling dates/times, including MJD, UTC, and TAI.
Definition: DateTime.h:58
#define PTR(...)
Definition: base.h:41
table::Key< std::string > name
Definition: ApCorrMap.cc:71
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
table::Key< double > fluxMag0Sigma
Definition: Calib.cc:487
void swap(ImageBase< PixelT > &a, ImageBase< PixelT > &b)
Definition: Image.cc:291
An object passed to Persistable::write to allow it to persist itself.
virtual void remove(std::string const &name)
Definition: PropertySet.cc:754
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:47
bool operator==(Calib const &rhs) const
Definition: Calib.cc:197
A custom container class for records, based on std::vector.
Definition: Catalog.h:94
static bool _throwOnNegativeFlux
Definition: Calib.h:146
#define CONST_PTR(...)
Definition: base.h:47
std::pair< double, double > getFluxMag0() const
Definition: Calib.cc:274
void setFluxMag0(double fluxMag0, double fluxMag0Sigma=0.0)
Definition: Calib.cc:257
void setExptime(double exptime)
Definition: Calib.cc:241
static void setThrowOnNegativeFlux(bool raiseException)
Definition: Calib.cc:138
double get(DateSystem system=MJD, Timescale scale=TAI) const
Definition: DateTime.cc:419
table::Key< boost::int64_t > midTime
Definition: Calib.cc:484
double getFlux(double const mag) const
Definition: Calib.cc:329
table::Key< double > expTime
Definition: Calib.cc:485
double getExptime() const
Definition: Calib.cc:249
boost::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Definition: Catalog.h:469
table::Key< table::Array< Kernel::Pixel > > image
Definition: FixedKernel.cc:117
def log
Definition: log.py:85
lsst::daf::base::DateTime getMidTime() const
Definition: Calib.cc:219
virtual void write(OutputArchiveHandle &handle) const
Write the object to one or more catalogs.
Definition: Calib.cc:537
size_type size() const
Return the size of the first dimension.
tbl::Schema schema
int stripCalibKeywords(boost::shared_ptr< lsst::daf::base::PropertySet > metadata)
Definition: Calib.cc:159
detail::SimpleInitializer< N > allocate(Vector< int, N > const &shape)
Create an expression that allocates uninitialized memory for an array.
#define LSST_EXCEPT(type,...)
Definition: Exception.h:46
A multidimensional strided array.
Definition: Array.h:47
Base class for all records.
Definition: BaseRecord.h:27
Iterator end() const
Return an Iterator to one past the end of the array.
Definition: ArrayBase.h:108
double getAsDouble(std::string const &name) const
Definition: PropertySet.cc:406
virtual std::string getPersistenceName() const
Return the unique name used to persist this object and look up its factory.
Definition: Calib.cc:535
Class for storing generic metadata.
Definition: PropertySet.h:82
std::string getAsString(std::string const &name) const
Definition: PropertySet.cc:444
afw::table::Key< double > b
table::Key< double > fluxMag0
Definition: Calib.cc:486
Interface for PropertySet class.
ExpressionTraits< Derived >::Iterator Iterator
Nested expression or element iterator.
double _fluxMag0Sigma
Definition: Calib.h:145
static bool getThrowOnNegativeFlux()
Definition: Calib.cc:148
lsst::daf::base::DateTime _midTime
Definition: Calib.h:142
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
bool exists(std::string const &name) const
Definition: PropertySet.cc:190
void setMidTime(lsst::daf::base::DateTime const &midTime)
Definition: Calib.cc:210
Include files required for standard LSST Exception handling.
Iterator begin() const
Return an Iterator to the beginning of the array.
Definition: ArrayBase.h:99