LSSTApplications  17.0+11,17.0+34,17.0+56,17.0+57,17.0+59,17.0+7,17.0-1-g377950a+33,17.0.1-1-g114240f+2,17.0.1-1-g4d4fbc4+28,17.0.1-1-g55520dc+49,17.0.1-1-g5f4ed7e+52,17.0.1-1-g6dd7d69+17,17.0.1-1-g8de6c91+11,17.0.1-1-gb9095d2+7,17.0.1-1-ge9fec5e+5,17.0.1-1-gf4e0155+55,17.0.1-1-gfc65f5f+50,17.0.1-1-gfc6fb1f+20,17.0.1-10-g87f9f3f+1,17.0.1-11-ge9de802+16,17.0.1-16-ga14f7d5c+4,17.0.1-17-gc79d625+1,17.0.1-17-gdae4c4a+8,17.0.1-2-g26618f5+29,17.0.1-2-g54f2ebc+9,17.0.1-2-gf403422+1,17.0.1-20-g2ca2f74+6,17.0.1-23-gf3eadeb7+1,17.0.1-3-g7e86b59+39,17.0.1-3-gb5ca14a,17.0.1-3-gd08d533+40,17.0.1-30-g596af8797,17.0.1-4-g59d126d+4,17.0.1-4-gc69c472+5,17.0.1-6-g5afd9b9+4,17.0.1-7-g35889ee+1,17.0.1-7-gc7c8782+18,17.0.1-9-gc4bbfb2+3,w.2019.22
LSSTDataManagementBasePackage
PhotoCalib.cc
Go to the documentation of this file.
1 /*
2  * LSST Data Management System
3  * Copyright 2017 LSST Corporation.
4  *
5  * This product includes software developed by the
6  * LSST Project (http://www.lsst.org/).
7  *
8  * This program is free software: you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation, either version 3 of the License, or
11  * (at your option) any later version.
12  *
13  * This program is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the LSST License Statement and
19  * the GNU General Public License along with this program. If not,
20  * see <http://www.lsstcorp.org/LegalNotices/>.
21  */
22 
23 #include <cmath>
24 #include <iostream>
25 
26 #include "lsst/geom/Point.h"
29 #include "lsst/afw/table.h"
33 #include "lsst/pex/exceptions.h"
34 #include "ndarray.h"
35 #include "lsst/afw/table/io/Persistable.cc"
36 #include "lsst/utils/Magnitude.h"
37 
38 namespace lsst {
39 namespace afw {
40 
43 
44 namespace image {
45 
46 // ------------------- helpers -------------------
47 
50  s << std::setprecision(std::numeric_limits<long double>::digits10 + 1);
51  s << "value=" << measurement.value << ", error=" << measurement.error;
52  return os << s.str();
53 }
54 
55 namespace {
56 
57 int const SERIALIZATION_VERSION = 1;
58 
59 double toNanojansky(double instFlux, double scale) { return instFlux * scale; }
60 
61 double toMagnitude(double instFlux, double scale) { return utils::nanojanskyToABMagnitude(instFlux * scale); }
62 
63 double toInstFluxFromMagnitude(double magnitude, double scale) {
64  // Note: flux[nJy] / scale = instFlux[counts]
65  return utils::ABMagnitudeToNanojansky(magnitude) / scale;
66 }
67 
68 double toNanojanskyErr(double instFlux, double instFluxErr, double scale, double scaleErr,
69  double nanojansky) {
70  return std::abs(nanojansky) * hypot(instFluxErr / instFlux, scaleErr / scale);
71 }
72 
86 void toNanojanskyVariance(ndarray::Array<float const, 2, 1> const &instFlux,
87  ndarray::Array<float const, 2, 1> const &instFluxVar, float scaleErr,
88  ndarray::Array<float const, 2, 1> const &flux, ndarray::Array<float, 2, 1> out) {
89  auto eigenFlux = ndarray::asEigen<Eigen::ArrayXpr>(flux);
90  auto eigenInstFluxVar = ndarray::asEigen<Eigen::ArrayXpr>(instFluxVar);
91  auto eigenInstFlux = ndarray::asEigen<Eigen::ArrayXpr>(instFlux);
92  auto eigenOut = ndarray::asEigen<Eigen::ArrayXpr>(out);
93  eigenOut = eigenFlux.square() *
94  (eigenInstFluxVar / eigenInstFlux.square() + (scaleErr / eigenFlux * eigenInstFlux).square());
95 }
96 
97 double toMagnitudeErr(double instFlux, double instFluxErr, double scale, double scaleErr) {
98  return 2.5 / log(10.0) * hypot(instFluxErr / instFlux, scaleErr / scale);
99 }
100 
101 } // anonymous namespace
102 
103 // ------------------ Deprecated Calib interface -------------------
104 
105 ndarray::Array<double, 1> PhotoCalib::getMagnitude(ndarray::Array<double const, 1> const &instFlux) const {
106  ndarray::Array<double, 1> result = ndarray::allocate(ndarray::makeVector(int(instFlux.size())));
107  auto iter = result.begin();
108  for (auto const &i : instFlux) {
109  *iter = toMagnitude(i, _calibrationMean);
110  iter++;
111  }
112  return result;
113 }
114 
115 std::pair<ndarray::Array<double, 1>, ndarray::Array<double, 1>> PhotoCalib::getMagnitude(
116  ndarray::Array<double const, 1> const &instFlux,
117  ndarray::Array<double const, 1> const &instFluxErr) const {
118  ndarray::Array<double, 1> mag = ndarray::allocate(ndarray::makeVector(int(instFlux.size())));
119  ndarray::Array<double, 1> magErr = ndarray::allocate(ndarray::makeVector(int(instFlux.size())));
120  auto iMag = mag.begin();
121  auto iMagErr = magErr.begin();
122  for (auto i = instFlux.begin(), iErr = instFluxErr.begin(); i != instFlux.end(); ++i, ++iErr) {
123  *iMag = toMagnitude(*i, _calibrationMean);
124  *iMagErr = toMagnitudeErr(*i, *iErr, _calibrationMean, _calibrationErr);
125  iMag++;
126  iMagErr++;
127  }
128  return std::make_pair(mag, magErr);
129 }
130 
131 // ------------------- Conversions to nanojansky -------------------
132 
133 double PhotoCalib::instFluxToNanojansky(double instFlux, lsst::geom::Point<double, 2> const &point) const {
134  return toNanojansky(instFlux, evaluate(point));
135 }
136 
137 double PhotoCalib::instFluxToNanojansky(double instFlux) const {
138  return toNanojansky(instFlux, _calibrationMean);
139 }
140 
141 Measurement PhotoCalib::instFluxToNanojansky(double instFlux, double instFluxErr,
142  lsst::geom::Point<double, 2> const &point) const {
143  double calibration, error, nanojansky;
144  calibration = evaluate(point);
145  nanojansky = toNanojansky(instFlux, calibration);
146  error = toNanojanskyErr(instFlux, instFluxErr, calibration, _calibrationErr, nanojansky);
147  return Measurement(nanojansky, error);
148 }
149 
150 Measurement PhotoCalib::instFluxToNanojansky(double instFlux, double instFluxErr) const {
151  double nanojansky = toNanojansky(instFlux, _calibrationMean);
152  double error = toNanojanskyErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr, nanojansky);
153  return Measurement(nanojansky, error);
154 }
155 
156 Measurement PhotoCalib::instFluxToNanojansky(afw::table::SourceRecord const &sourceRecord,
157  std::string const &instFluxField) const {
158  auto position = sourceRecord.getCentroid();
159  auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
160  auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
161  return instFluxToNanojansky(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
162 }
163 ndarray::Array<double, 2, 2> PhotoCalib::instFluxToNanojansky(afw::table::SourceCatalog const &sourceCatalog,
164  std::string const &instFluxField) const {
165  ndarray::Array<double, 2, 2> result =
166  ndarray::allocate(ndarray::makeVector(int(sourceCatalog.size()), 2));
167  instFluxToNanojanskyArray(sourceCatalog, instFluxField, result);
168  return result;
169 }
170 
171 void PhotoCalib::instFluxToNanojansky(afw::table::SourceCatalog &sourceCatalog,
172  std::string const &instFluxField, std::string const &outField) const {
173  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
174  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
175  auto nanojanskyKey = sourceCatalog.getSchema().find<double>(outField + "_flux").key;
176  auto nanojanskyErrKey = sourceCatalog.getSchema().find<double>(outField + "_fluxErr").key;
177  for (auto &record : sourceCatalog) {
178  auto result = instFluxToNanojansky(record.get(instFluxKey), record.get(instFluxErrKey),
179  record.getCentroid());
180  record.set(nanojanskyKey, result.value);
181  record.set(nanojanskyErrKey, result.error);
182  }
183 }
184 
185 // ------------------- Conversions to Magnitudes -------------------
186 
187 double PhotoCalib::instFluxToMagnitude(double instFlux, lsst::geom::Point<double, 2> const &point) const {
188  return toMagnitude(instFlux, evaluate(point));
189 }
190 
191 double PhotoCalib::instFluxToMagnitude(double instFlux) const {
192  return toMagnitude(instFlux, _calibrationMean);
193 }
194 
195 Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr,
196  lsst::geom::Point<double, 2> const &point) const {
197  double calibration, error, magnitude;
198  calibration = evaluate(point);
199  magnitude = toMagnitude(instFlux, calibration);
200  error = toMagnitudeErr(instFlux, instFluxErr, calibration, _calibrationErr);
201  return Measurement(magnitude, error);
202 }
203 
204 Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr) const {
205  double magnitude = toMagnitude(instFlux, _calibrationMean);
206  double error = toMagnitudeErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr);
207  return Measurement(magnitude, error);
208 }
209 
210 Measurement PhotoCalib::instFluxToMagnitude(afw::table::SourceRecord const &sourceRecord,
211  std::string const &instFluxField) const {
212  auto position = sourceRecord.getCentroid();
213  auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
214  auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
215  return instFluxToMagnitude(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
216 }
217 
218 ndarray::Array<double, 2, 2> PhotoCalib::instFluxToMagnitude(afw::table::SourceCatalog const &sourceCatalog,
219  std::string const &instFluxField) const {
220  ndarray::Array<double, 2, 2> result =
221  ndarray::allocate(ndarray::makeVector(int(sourceCatalog.size()), 2));
222  instFluxToMagnitudeArray(sourceCatalog, instFluxField, result);
223  return result;
224 }
225 
226 void PhotoCalib::instFluxToMagnitude(afw::table::SourceCatalog &sourceCatalog,
227  std::string const &instFluxField, std::string const &outField) const {
228  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
229  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
230  auto magKey = sourceCatalog.getSchema().find<double>(outField + "_mag").key;
231  auto magErrKey = sourceCatalog.getSchema().find<double>(outField + "_magErr").key;
232  for (auto &record : sourceCatalog) {
233  auto result = instFluxToMagnitude(record.get(instFluxKey), record.get(instFluxErrKey),
234  record.getCentroid());
235  record.set(magKey, result.value);
236  record.set(magErrKey, result.error);
237  }
238 }
239 
240 // ------------------- other utility methods -------------------
241 
242 double PhotoCalib::magnitudeToInstFlux(double magnitude) const {
243  return toInstFluxFromMagnitude(magnitude, _calibrationMean);
244 }
245 
246 double PhotoCalib::magnitudeToInstFlux(double magnitude, lsst::geom::Point<double, 2> const &point) const {
247  return toInstFluxFromMagnitude(magnitude, evaluate(point));
248 }
249 
250 std::shared_ptr<math::BoundedField> PhotoCalib::computeScaledCalibration() const {
251  return *(_calibration) / _calibrationMean;
252 }
253 
255  throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented: See DM-10154.");
256 }
257 
258 bool PhotoCalib::operator==(PhotoCalib const &rhs) const {
259  return (_calibrationMean == rhs._calibrationMean && _calibrationErr == rhs._calibrationErr &&
260  (*_calibration) == *(rhs._calibration));
261 }
262 
263 double PhotoCalib::computeCalibrationMean(std::shared_ptr<afw::math::BoundedField> calibration) const {
264  return calibration->mean();
265 }
266 
267 std::shared_ptr<typehandling::Storable> PhotoCalib::cloneStorable() const {
268  return std::make_unique<PhotoCalib>(*this);
269 }
270 
271 std::string PhotoCalib::toString() const {
272  std::stringstream buffer;
273  if (_isConstant)
274  buffer << "spatially constant with ";
275  else
276  buffer << *_calibration << " with ";
277  buffer << "mean: " << _calibrationMean << " error: " << _calibrationErr;
278  return buffer.str();
279 }
280 
281 bool PhotoCalib::equals(typehandling::Storable const &other) const noexcept {
282  return singleClassEquals(*this, other);
283 }
284 
286  return os << photoCalib.toString();
287 }
288 
289 MaskedImage<float> PhotoCalib::calibrateImage(MaskedImage<float> const &maskedImage,
290  bool includeScaleUncertainty) const {
291  // Deep copy construct, as we're mutiplying in-place.
292  auto result = MaskedImage<float>(maskedImage, true);
293 
294  if (_isConstant) {
295  *(result.getImage()) *= _calibrationMean;
296  } else {
297  _calibration->multiplyImage(*(result.getImage()), true); // only in the overlap region
298  }
299  if (includeScaleUncertainty) {
300  toNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(),
301  _calibrationErr, result.getImage()->getArray(),
302  result.getVariance()->getArray());
303  } else {
304  toNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(), 0,
305  result.getImage()->getArray(), result.getVariance()->getArray());
306  }
307 
308  return result;
309 }
310 
311 afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog,
312  std::vector<std::string> const &instFluxFields) const {
313  auto const &inSchema = catalog.getSchema();
314  afw::table::SchemaMapper mapper(inSchema, true); // true: share the alias map
315  mapper.addMinimalSchema(inSchema);
316 
317  using FieldD = afw::table::Field<double>;
318 
319  struct Keys {
320  table::Key<double> instFlux;
321  table::Key<double> instFluxErr;
322  table::Key<double> flux;
323  table::Key<double> fluxErr;
324  table::Key<double> mag;
325  table::Key<double> magErr;
326  };
327 
329  keys.reserve(instFluxFields.size());
330  for (auto const &field : instFluxFields) {
331  Keys newKey;
332  newKey.instFlux = inSchema[inSchema.join(field, "instFlux")];
333  newKey.flux =
334  mapper.addOutputField(FieldD(inSchema.join(field, "flux"), "calibrated flux", "nJy"), true);
335  newKey.mag = mapper.addOutputField(
336  FieldD(inSchema.join(field, "mag"), "calibrated magnitude", "mag(AB)"), true);
337  try {
338  newKey.instFluxErr = inSchema.find<double>(inSchema.join(field, "instFluxErr")).key;
339  newKey.fluxErr = mapper.addOutputField(
340  FieldD(inSchema.join(field, "fluxErr"), "calibrated flux uncertainty", "nJy"), true);
341  newKey.magErr = mapper.addOutputField(
342  FieldD(inSchema.join(field, "magErr"), "calibrated magnitude uncertainty", "mag(AB)"),
343  true);
344  } catch (pex::exceptions::NotFoundError &) {
345  ; // Keys struct defaults to invalid keys; that marks the error as missing.
346  }
347  keys.emplace_back(newKey);
348  }
349 
350  // Create the new catalog
352  output.insert(mapper, output.begin(), catalog.begin(), catalog.end());
353 
354  auto calibration = evaluateCatalog(output);
355 
356  // fill in the catalog values
357  int iRec = 0;
358  for (auto &rec : output) {
359  for (auto &key : keys) {
360  double instFlux = rec.get(key.instFlux);
361  double nanojansky = toNanojansky(instFlux, calibration[iRec]);
362  rec.set(key.flux, nanojansky);
363  rec.set(key.mag, toMagnitude(instFlux, calibration[iRec]));
364  if (key.instFluxErr.isValid()) {
365  double instFluxErr = rec.get(key.instFluxErr);
366  rec.set(key.fluxErr, toNanojanskyErr(instFlux, instFluxErr, calibration[iRec],
367  _calibrationErr, nanojansky));
368  rec.set(key.magErr,
369  toMagnitudeErr(instFlux, instFluxErr, calibration[iRec], _calibrationErr));
370  }
371  }
372  ++iRec;
373  }
374 
375  return output;
376 }
377 
378 afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog) const {
379  std::vector<std::string> instFluxFields;
380  static std::string const SUFFIX = "_instFlux";
381  for (auto const &name : catalog.getSchema().getNames()) {
382  // Pick every field ending in "_instFlux", grabbing everything before that prefix.
383  if (name.size() > SUFFIX.size() + 1 &&
384  name.compare(name.size() - SUFFIX.size(), SUFFIX.size(), SUFFIX) == 0) {
385  instFluxFields.emplace_back(name.substr(0, name.size() - 9));
386  }
387  }
388  return calibrateCatalog(catalog, instFluxFields);
389 }
390 
391 // ------------------- persistence -------------------
392 
393 namespace {
394 
395 class PhotoCalibSchema {
396 public:
397  table::Schema schema;
398  table::Key<double> calibrationMean;
399  table::Key<double> calibrationErr;
400  table::Key<table::Flag> isConstant;
401  table::Key<int> field;
402  table::Key<int> version;
403 
404  // No copying
405  PhotoCalibSchema(PhotoCalibSchema const &) = delete;
406  PhotoCalibSchema &operator=(PhotoCalibSchema const &) = delete;
407  // No moving
408  PhotoCalibSchema(PhotoCalibSchema &&) = delete;
409  PhotoCalibSchema &operator=(PhotoCalibSchema &&) = delete;
410 
411  static PhotoCalibSchema const &get() {
412  static PhotoCalibSchema const instance;
413  return instance;
414  }
415 
416 private:
417  PhotoCalibSchema()
418  : schema(),
419  calibrationMean(schema.addField<double>(
420  "calibrationMean", "mean calibration on this PhotoCalib's domain", "count")),
422  schema.addField<double>("calibrationErr", "1-sigma error on calibrationMean", "count")),
423  isConstant(schema.addField<table::Flag>("isConstant", "Is this spatially-constant?")),
424  field(schema.addField<int>("field", "archive ID of the BoundedField object")),
425  version(schema.addField<int>("version", "version of this PhotoCalib")) {
426  schema.getCitizen().markPersistent();
427  }
428 };
429 
430 class PhotoCalibFactory : public table::io::PersistableFactory {
431 public:
432  PTR(table::io::Persistable)
433  read(InputArchive const &archive, CatalogVector const &catalogs) const override {
434  table::BaseRecord const &record = catalogs.front().front();
435  PhotoCalibSchema const &keys = PhotoCalibSchema::get();
436  int version = getVersion(record);
437  if (version < 1) {
438  throw(pex::exceptions::RuntimeError("Unsupported version (version 0 was defined in maggies): " +
439  std::to_string(version)));
440  }
441  return std::make_shared<PhotoCalib>(record.get(keys.calibrationMean), record.get(keys.calibrationErr),
442  archive.get<afw::math::BoundedField>(record.get(keys.field)),
443  record.get(keys.isConstant));
444  }
445 
446  PhotoCalibFactory(std::string const &name) : afw::table::io::PersistableFactory(name) {}
447 
448 protected:
449  int getVersion(table::BaseRecord const &record) const {
450  int version = -1;
451  try {
452  std::string versionName = "version";
453  auto versionKey = record.getSchema().find<int>(versionName);
454  version = record.get(versionKey.key);
455  } catch (const pex::exceptions::NotFoundError &) {
456  // un-versioned files are version 0
457  version = 0;
458  }
459  return version;
460  }
461 };
462 
463 std::string getPhotoCalibPersistenceName() { return "PhotoCalib"; }
464 
465 PhotoCalibFactory registration(getPhotoCalibPersistenceName());
466 
467 } // namespace
468 
473 namespace {
474 int const CALIB_TABLE_CURRENT_VERSION = 2; // final version of Calib in ExposureTable
475 std::string const EXPTIME_FIELD_NAME = "exptime"; // name of exposure time field (no longer used)
476 
477 class CalibKeys {
478 public:
479  table::Schema schema;
480  table::Key<std::int64_t> midTime;
481  table::Key<double> expTime;
482  table::Key<double> fluxMag0;
483  table::Key<double> fluxMag0Err;
484 
485  // No copying
486  CalibKeys(const CalibKeys &) = delete;
487  CalibKeys &operator=(const CalibKeys &) = delete;
488 
489  // No moving
490  CalibKeys(CalibKeys &&) = delete;
491  CalibKeys &operator=(CalibKeys &&) = delete;
492 
493  CalibKeys(int tableVersion = CALIB_TABLE_CURRENT_VERSION)
494  : schema(), midTime(), expTime(), fluxMag0(), fluxMag0Err() {
495  if (tableVersion == 1) {
496  // obsolete fields
497  midTime = schema.addField<std::int64_t>(
498  "midtime", "middle of the time of the exposure relative to Unix epoch", "ns");
499  expTime = schema.addField<double>(EXPTIME_FIELD_NAME, "exposure time", "s");
500  }
501  fluxMag0 = schema.addField<double>("fluxmag0", "flux of a zero-magnitude object", "count");
502  fluxMag0Err = schema.addField<double>("fluxmag0.err", "1-sigma error on fluxmag0", "count");
503  }
504 };
505 
506 class CalibFactory : public table::io::PersistableFactory {
507 public:
508  std::shared_ptr<table::io::Persistable> read(InputArchive const &archive,
509  CatalogVector const &catalogs) const override {
510  // table version is not persisted, so we don't have a clean way to determine the version;
511  // the hack is version = 1 if exptime found, else current
512  int tableVersion = 1;
513  try {
514  catalogs.front().getSchema().find<double>(EXPTIME_FIELD_NAME);
515  } catch (pex::exceptions::NotFoundError) {
516  tableVersion = CALIB_TABLE_CURRENT_VERSION;
517  }
518 
519  CalibKeys const keys{tableVersion};
520  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
521  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
522  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
523  table::BaseRecord const &record = catalogs.front().front();
524 
525  double calibration = utils::referenceFlux / record.get(keys.fluxMag0);
526  double calibrationErr =
527  utils::referenceFlux * record.get(keys.fluxMag0Err) / std::pow(record.get(keys.fluxMag0), 2);
528  return std::make_shared<PhotoCalib>(calibration, calibrationErr);
529  }
530 
531  explicit CalibFactory(std::string const &name) : table::io::PersistableFactory(name) {}
532 };
533 
534 std::string getCalibPersistenceName() { return "Calib"; }
535 
536 CalibFactory calibRegistration(getCalibPersistenceName());
537 
538 } // namespace
539 
540 std::string PhotoCalib::getPersistenceName() const { return getPhotoCalibPersistenceName(); }
541 
542 void PhotoCalib::write(OutputArchiveHandle &handle) const {
543  PhotoCalibSchema const &keys = PhotoCalibSchema::get();
544  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
545  auto record = catalog.addNew();
546  record->set(keys.calibrationMean, _calibrationMean);
547  record->set(keys.calibrationErr, _calibrationErr);
548  record->set(keys.isConstant, _isConstant);
549  record->set(keys.field, handle.put(_calibration));
550  record->set(keys.version, SERIALIZATION_VERSION);
551  handle.saveCatalog(catalog);
552 }
553 
554 // ------------------- private/protected helpers -------------------
555 
556 double PhotoCalib::evaluate(lsst::geom::Point<double, 2> const &point) const {
557  if (_isConstant)
558  return _calibrationMean;
559  else
560  return _calibration->evaluate(point);
561 }
562 
563 ndarray::Array<double, 1> PhotoCalib::evaluateArray(ndarray::Array<double, 1> const &xx,
564  ndarray::Array<double, 1> const &yy) const {
565  if (_isConstant) {
566  ndarray::Array<double, 1> result = ndarray::allocate(ndarray::makeVector(xx.size()));
567  result.deep() = _calibrationMean;
568  return result;
569  } else {
570  return _calibration->evaluate(xx, yy);
571  }
572 }
573 
574 ndarray::Array<double, 1> PhotoCalib::evaluateCatalog(afw::table::SourceCatalog const &sourceCatalog) const {
575  ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
576  ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
577  size_t i = 0;
578  for (auto const &rec : sourceCatalog) {
579  auto point = rec.getCentroid();
580  xx[i] = point.getX();
581  yy[i] = point.getY();
582  ++i;
583  }
584  return evaluateArray(xx, yy);
585 }
586 
587 void PhotoCalib::instFluxToNanojanskyArray(afw::table::SourceCatalog const &sourceCatalog,
588  std::string const &instFluxField,
589  ndarray::Array<double, 2, 2> result) const {
590  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
591  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
592 
593  auto calibration = evaluateCatalog(sourceCatalog);
594  int i = 0;
595  auto iter = result.begin();
596  for (auto const &rec : sourceCatalog) {
597  double instFlux = rec.get(instFluxKey);
598  double instFluxErr = rec.get(instFluxErrKey);
599  double nanojansky = toNanojansky(instFlux, calibration[i]);
600  (*iter)[0] = nanojansky;
601  (*iter)[1] = toNanojanskyErr(instFlux, instFluxErr, calibration[i], _calibrationErr, nanojansky);
602  ++iter;
603  ++i;
604  }
605 }
606 
607 void PhotoCalib::instFluxToMagnitudeArray(afw::table::SourceCatalog const &sourceCatalog,
608  std::string const &instFluxField,
609  ndarray::Array<double, 2, 2> result) const {
610  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
611  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
612 
613  auto calibration = evaluateCatalog(sourceCatalog);
614  auto iter = result.begin();
615  int i = 0;
616  for (auto const &rec : sourceCatalog) {
617  double instFlux = rec.get(instFluxKey);
618  double instFluxErr = rec.get(instFluxErrKey);
619  (*iter)[0] = toMagnitude(instFlux, calibration[i]);
620  (*iter)[1] = toMagnitudeErr(instFlux, instFluxErr, calibration[i], _calibrationErr);
621  ++iter;
622  ++i;
623  }
624 }
625 
627  auto key = "FLUXMAG0";
628  if (metadata.exists(key)) {
629  double instFluxMag0 = metadata.getAsDouble(key);
630  if (strip) metadata.remove(key);
631 
632  double instFluxMag0Err = 0.0;
633  key = "FLUXMAG0ERR";
634  if (metadata.exists(key)) {
635  instFluxMag0Err = metadata.getAsDouble(key);
636  if (strip) metadata.remove(key);
637  }
638  return makePhotoCalibFromCalibZeroPoint(instFluxMag0, instFluxMag0Err);
639  } else {
640  return nullptr;
641  }
642 }
643 
644 std::shared_ptr<PhotoCalib> makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err) {
645  double calibration = utils::referenceFlux / instFluxMag0;
646  double calibrationErr = utils::referenceFlux * instFluxMag0Err / std::pow(instFluxMag0, 2);
647  return std::make_shared<PhotoCalib>(calibration, calibrationErr);
648 }
649 
650 } // namespace image
651 } // namespace afw
652 } // namespace lsst
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage&#39;s variance.
Definition: MaskedImage.h:1091
std::shared_ptr< PhotoCalib > makePhotoCalibFromMetadata(daf::base::PropertySet &metadata, bool strip=false)
Construct a PhotoCalib from FITS FLUXMAG0/FLUXMAG0ERR keywords.
Definition: PhotoCalib.cc:626
std::string toString() const override
Create a string representation of this object.
Definition: PhotoCalib.cc:271
Angle abs(Angle const &a)
Definition: Angle.h:106
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
Definition: Persistable.cc:18
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values...
Definition: PhotoCalib.cc:644
int put(Persistable const *obj, bool permissive=false)
Save an object to the archive and return a unique ID that can be used to retrieve it from an InputArc...
table::Key< double > fluxMag0
Definition: PhotoCalib.cc:482
An object passed to Persistable::write to allow it to persist itself.
The photometric calibration of an exposure.
Definition: PhotoCalib.h:116
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
Interface supporting iteration over heterogenous containers.
Definition: Storable.h:56
Key< T > addOutputField(Field< T > const &newField, bool doReplace=false)
Add a new field to the output Schema that is not connected to the input Schema.
Definition: SchemaMapper.h:34
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:109
T to_string(T... args)
Schema const getOutputSchema() const
Return the output schema (copy-on-write).
Definition: SchemaMapper.h:27
py::object result
Definition: schema.cc:418
table::Key< int > field
Definition: PhotoCalib.cc:401
A value and its error.
Definition: PhotoCalib.h:53
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage&#39;s image.
Definition: MaskedImage.h:1058
STL class.
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition: Schema.cc:656
A base class for image defects.
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:74
Utilities for converting between flux and magnitude in C++.
A description of a field in a table.
Definition: Field.h:24
T str(T... args)
T make_pair(T... args)
Schema getSchema() const
Return the Schema that holds this record&#39;s fields and keys.
Definition: BaseRecord.h:80
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
table::Key< double > calibrationMean
Definition: PhotoCalib.cc:398
table::Key< double > expTime
Definition: PhotoCalib.cc:481
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
solver_t * s
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
std::ostream & operator<<(std::ostream &os, PhotoCalib const &photoCalib)
Definition: PhotoCalib.cc:285
T find(T... args)
T size(T... args)
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
double nanojanskyToABMagnitude(double flux)
Convert a flux in nanojansky to AB magnitude.
Definition: Magnitude.cc:30
table::Schema schema
Definition: PhotoCalib.cc:397
void addMinimalSchema(Schema const &minimal, bool doMap=true)
Add the given minimal schema to the output schema.
Key< U > key
Definition: Schema.cc:281
T pow(T... args)
table::Key< std::int64_t > midTime
Definition: PhotoCalib.cc:480
Class for storing generic metadata.
Definition: PropertySet.h:68
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition: Persistable.h:48
#define PTR(...)
Definition: base.h:41
table::Key< double > fluxMag0Err
Definition: PhotoCalib.cc:483
ItemVariant const * other
Definition: Schema.cc:56
Record class that contains measurements made on a single exposure.
Definition: Source.h:82
table::Key< table::Flag > isConstant
Definition: PhotoCalib.cc:400
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
table::Key< int > version
Definition: PhotoCalib.cc:402
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects...
const double referenceFlux
The Oke & Gunn (1983) AB magnitude reference flux, in nJy (often approximated as 3631.0).
Definition: Magnitude.h:46
STL class.
double ABMagnitudeToNanojansky(double magnitude)
Convert an AB magnitude to a flux in nanojansky.
Definition: Magnitude.cc:32
Implementation of the Photometric Calibration class.
bool strip
Definition: fits.cc:883
table::Key< double > calibrationErr
Definition: PhotoCalib.cc:399
std::ostream * os
Definition: Schema.cc:746
SchemaMapper * mapper
Definition: SchemaMapper.cc:78
T reserve(T... args)
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
Definition: Source.h:656
std::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:485
T emplace_back(T... args)