LSST Applications  21.0.0+04719a4bac,21.0.0-1-ga51b5d4+f5e6047307,21.0.0-11-g2b59f77+a9c1acf22d,21.0.0-11-ga42c5b2+86977b0b17,21.0.0-12-gf4ce030+76814010d2,21.0.0-13-g1721dae+760e7a6536,21.0.0-13-g3a573fe+768d78a30a,21.0.0-15-g5a7caf0+f21cbc5713,21.0.0-16-g0fb55c1+b60e2d390c,21.0.0-19-g4cded4ca+71a93a33c0,21.0.0-2-g103fe59+bb20972958,21.0.0-2-g45278ab+04719a4bac,21.0.0-2-g5242d73+3ad5d60fb1,21.0.0-2-g7f82c8f+8babb168e8,21.0.0-2-g8f08a60+06509c8b61,21.0.0-2-g8faa9b5+616205b9df,21.0.0-2-ga326454+8babb168e8,21.0.0-2-gde069b7+5e4aea9c2f,21.0.0-2-gecfae73+1d3a86e577,21.0.0-2-gfc62afb+3ad5d60fb1,21.0.0-25-g1d57be3cd+e73869a214,21.0.0-3-g357aad2+ed88757d29,21.0.0-3-g4a4ce7f+3ad5d60fb1,21.0.0-3-g4be5c26+3ad5d60fb1,21.0.0-3-g65f322c+e0b24896a3,21.0.0-3-g7d9da8d+616205b9df,21.0.0-3-ge02ed75+a9c1acf22d,21.0.0-4-g591bb35+a9c1acf22d,21.0.0-4-g65b4814+b60e2d390c,21.0.0-4-gccdca77+0de219a2bc,21.0.0-4-ge8a399c+6c55c39e83,21.0.0-5-gd00fb1e+05fce91b99,21.0.0-6-gc675373+3ad5d60fb1,21.0.0-64-g1122c245+4fb2b8f86e,21.0.0-7-g04766d7+cd19d05db2,21.0.0-7-gdf92d54+04719a4bac,21.0.0-8-g5674e7b+d1bd76f71f,master-gac4afde19b+a9c1acf22d,w.2021.13
LSST Data Management Base Package
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"
36 #include "lsst/utils/Magnitude.h"
37 
38 namespace lsst {
39 namespace afw {
40 
43 
44 namespace image {
45 
46 // ------------------- helpers -------------------
47 
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 // ------------------- Conversions to nanojansky -------------------
104 
105 double PhotoCalib::instFluxToNanojansky(double instFlux, lsst::geom::Point<double, 2> const &point) const {
106  return toNanojansky(instFlux, evaluate(point));
107 }
108 
109 double PhotoCalib::instFluxToNanojansky(double instFlux) const {
110  return toNanojansky(instFlux, _calibrationMean);
111 }
112 
113 Measurement PhotoCalib::instFluxToNanojansky(double instFlux, double instFluxErr,
114  lsst::geom::Point<double, 2> const &point) const {
115  double calibration, error, nanojansky;
116  calibration = evaluate(point);
117  nanojansky = toNanojansky(instFlux, calibration);
118  error = toNanojanskyErr(instFlux, instFluxErr, calibration, _calibrationErr, nanojansky);
119  return Measurement(nanojansky, error);
120 }
121 
122 Measurement PhotoCalib::instFluxToNanojansky(double instFlux, double instFluxErr) const {
123  double nanojansky = toNanojansky(instFlux, _calibrationMean);
124  double error = toNanojanskyErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr, nanojansky);
125  return Measurement(nanojansky, error);
126 }
127 
128 Measurement PhotoCalib::instFluxToNanojansky(afw::table::SourceRecord const &sourceRecord,
129  std::string const &instFluxField) const {
130  auto position = sourceRecord.getCentroid();
131  auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
132  auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
133  return instFluxToNanojansky(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
134 }
135 ndarray::Array<double, 2, 2> PhotoCalib::instFluxToNanojansky(afw::table::SourceCatalog const &sourceCatalog,
136  std::string const &instFluxField) const {
137  ndarray::Array<double, 2, 2> result =
138  ndarray::allocate(ndarray::makeVector(int(sourceCatalog.size()), 2));
139  instFluxToNanojanskyArray(sourceCatalog, instFluxField, result);
140  return result;
141 }
142 
143 void PhotoCalib::instFluxToNanojansky(afw::table::SourceCatalog &sourceCatalog,
144  std::string const &instFluxField, std::string const &outField) const {
145  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
146  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
147  auto nanojanskyKey = sourceCatalog.getSchema().find<double>(outField + "_flux").key;
148  auto nanojanskyErrKey = sourceCatalog.getSchema().find<double>(outField + "_fluxErr").key;
149  for (auto &record : sourceCatalog) {
150  auto result = instFluxToNanojansky(record.get(instFluxKey), record.get(instFluxErrKey),
151  record.getCentroid());
152  record.set(nanojanskyKey, result.value);
153  record.set(nanojanskyErrKey, result.error);
154  }
155 }
156 
157 // ------------------- Conversions to Magnitudes -------------------
158 
159 double PhotoCalib::instFluxToMagnitude(double instFlux, lsst::geom::Point<double, 2> const &point) const {
160  return toMagnitude(instFlux, evaluate(point));
161 }
162 
163 double PhotoCalib::instFluxToMagnitude(double instFlux) const {
164  return toMagnitude(instFlux, _calibrationMean);
165 }
166 
167 Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr,
168  lsst::geom::Point<double, 2> const &point) const {
169  double calibration, error, magnitude;
170  calibration = evaluate(point);
171  magnitude = toMagnitude(instFlux, calibration);
172  error = toMagnitudeErr(instFlux, instFluxErr, calibration, _calibrationErr);
173  return Measurement(magnitude, error);
174 }
175 
176 Measurement PhotoCalib::instFluxToMagnitude(double instFlux, double instFluxErr) const {
177  double magnitude = toMagnitude(instFlux, _calibrationMean);
178  double error = toMagnitudeErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr);
179  return Measurement(magnitude, error);
180 }
181 
182 Measurement PhotoCalib::instFluxToMagnitude(afw::table::SourceRecord const &sourceRecord,
183  std::string const &instFluxField) const {
184  auto position = sourceRecord.getCentroid();
185  auto instFluxKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFlux").key;
186  auto instFluxErrKey = sourceRecord.getSchema().find<double>(instFluxField + "_instFluxErr").key;
187  return instFluxToMagnitude(sourceRecord.get(instFluxKey), sourceRecord.get(instFluxErrKey), position);
188 }
189 
190 ndarray::Array<double, 2, 2> PhotoCalib::instFluxToMagnitude(afw::table::SourceCatalog const &sourceCatalog,
191  std::string const &instFluxField) const {
192  ndarray::Array<double, 2, 2> result =
193  ndarray::allocate(ndarray::makeVector(int(sourceCatalog.size()), 2));
194  instFluxToMagnitudeArray(sourceCatalog, instFluxField, result);
195  return result;
196 }
197 
198 void PhotoCalib::instFluxToMagnitude(afw::table::SourceCatalog &sourceCatalog,
199  std::string const &instFluxField, std::string const &outField) const {
200  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
201  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
202  auto magKey = sourceCatalog.getSchema().find<double>(outField + "_mag").key;
203  auto magErrKey = sourceCatalog.getSchema().find<double>(outField + "_magErr").key;
204  for (auto &record : sourceCatalog) {
205  auto result = instFluxToMagnitude(record.get(instFluxKey), record.get(instFluxErrKey),
206  record.getCentroid());
207  record.set(magKey, result.value);
208  record.set(magErrKey, result.error);
209  }
210 }
211 
212 // ------------------- other utility methods -------------------
213 
214 double PhotoCalib::magnitudeToInstFlux(double magnitude) const {
215  return toInstFluxFromMagnitude(magnitude, _calibrationMean);
216 }
217 
218 double PhotoCalib::magnitudeToInstFlux(double magnitude, lsst::geom::Point<double, 2> const &point) const {
219  return toInstFluxFromMagnitude(magnitude, evaluate(point));
220 }
221 
222 std::shared_ptr<math::BoundedField> PhotoCalib::computeScaledCalibration() const {
223  return *(_calibration) / _calibrationMean;
224 }
225 
227  throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented: See DM-10154.");
228 }
229 
230 bool PhotoCalib::operator==(PhotoCalib const &rhs) const {
231  return (_calibrationMean == rhs._calibrationMean && _calibrationErr == rhs._calibrationErr &&
232  (*_calibration) == *(rhs._calibration));
233 }
234 
235 double PhotoCalib::computeCalibrationMean(std::shared_ptr<afw::math::BoundedField> calibration) const {
236  return calibration->mean();
237 }
238 
239 std::shared_ptr<typehandling::Storable> PhotoCalib::cloneStorable() const {
240  return std::make_unique<PhotoCalib>(*this);
241 }
242 
243 std::string PhotoCalib::toString() const {
244  std::stringstream buffer;
245  if (_isConstant)
246  buffer << "spatially constant with ";
247  else
248  buffer << *_calibration << " with ";
249  buffer << "mean: " << _calibrationMean << " error: " << _calibrationErr;
250  return buffer.str();
251 }
252 
253 bool PhotoCalib::equals(typehandling::Storable const &other) const noexcept {
254  return singleClassEquals(*this, other);
255 }
256 
258  return os << photoCalib.toString();
259 }
260 
261 MaskedImage<float> PhotoCalib::calibrateImage(MaskedImage<float> const &maskedImage,
262  bool includeScaleUncertainty) const {
263  // Deep copy construct, as we're mutiplying in-place.
264  auto result = MaskedImage<float>(maskedImage, true);
265 
266  if (_isConstant) {
267  *(result.getImage()) *= _calibrationMean;
268  } else {
269  _calibration->multiplyImage(*(result.getImage()), true); // only in the overlap region
270  }
271  if (includeScaleUncertainty) {
272  toNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(),
273  _calibrationErr, result.getImage()->getArray(),
274  result.getVariance()->getArray());
275  } else {
276  toNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(), 0,
277  result.getImage()->getArray(), result.getVariance()->getArray());
278  }
279 
280  return result;
281 }
282 
283 afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog,
284  std::vector<std::string> const &instFluxFields) const {
285  auto const &inSchema = catalog.getSchema();
286  afw::table::SchemaMapper mapper(inSchema, true); // true: share the alias map
287  mapper.addMinimalSchema(inSchema);
288 
289  using FieldD = afw::table::Field<double>;
290 
291  struct Keys {
292  table::Key<double> instFlux;
293  table::Key<double> instFluxErr;
295  table::Key<double> fluxErr;
296  table::Key<double> mag;
297  table::Key<double> magErr;
298  };
299 
301  keys.reserve(instFluxFields.size());
302  for (auto const &field : instFluxFields) {
303  Keys newKey;
304  newKey.instFlux = inSchema[inSchema.join(field, "instFlux")];
305  newKey.flux =
306  mapper.addOutputField(FieldD(inSchema.join(field, "flux"), "calibrated flux", "nJy"), true);
307  newKey.mag = mapper.addOutputField(
308  FieldD(inSchema.join(field, "mag"), "calibrated magnitude", "mag(AB)"), true);
309  try {
310  newKey.instFluxErr = inSchema.find<double>(inSchema.join(field, "instFluxErr")).key;
311  newKey.fluxErr = mapper.addOutputField(
312  FieldD(inSchema.join(field, "fluxErr"), "calibrated flux uncertainty", "nJy"), true);
313  newKey.magErr = mapper.addOutputField(
314  FieldD(inSchema.join(field, "magErr"), "calibrated magnitude uncertainty", "mag(AB)"),
315  true);
316  } catch (pex::exceptions::NotFoundError &) {
317  ; // Keys struct defaults to invalid keys; that marks the error as missing.
318  }
319  keys.emplace_back(newKey);
320  }
321 
322  // Create the new catalog
323  afw::table::SourceCatalog output(mapper.getOutputSchema());
324  output.insert(mapper, output.begin(), catalog.begin(), catalog.end());
325 
326  auto calibration = evaluateCatalog(output);
327 
328  // fill in the catalog values
329  int iRec = 0;
330  for (auto &rec : output) {
331  for (auto &key : keys) {
332  double instFlux = rec.get(key.instFlux);
333  double nanojansky = toNanojansky(instFlux, calibration[iRec]);
334  rec.set(key.flux, nanojansky);
335  rec.set(key.mag, toMagnitude(instFlux, calibration[iRec]));
336  if (key.instFluxErr.isValid()) {
337  double instFluxErr = rec.get(key.instFluxErr);
338  rec.set(key.fluxErr, toNanojanskyErr(instFlux, instFluxErr, calibration[iRec],
339  _calibrationErr, nanojansky));
340  rec.set(key.magErr,
341  toMagnitudeErr(instFlux, instFluxErr, calibration[iRec], _calibrationErr));
342  }
343  }
344  ++iRec;
345  }
346 
347  return output;
348 }
349 
350 afw::table::SourceCatalog PhotoCalib::calibrateCatalog(afw::table::SourceCatalog const &catalog) const {
351  std::vector<std::string> instFluxFields;
352  static std::string const SUFFIX = "_instFlux";
353  for (auto const &name : catalog.getSchema().getNames()) {
354  // Pick every field ending in "_instFlux", grabbing everything before that prefix.
355  if (name.size() > SUFFIX.size() + 1 &&
356  name.compare(name.size() - SUFFIX.size(), SUFFIX.size(), SUFFIX) == 0) {
357  instFluxFields.emplace_back(name.substr(0, name.size() - 9));
358  }
359  }
360  return calibrateCatalog(catalog, instFluxFields);
361 }
362 
363 // ------------------- persistence -------------------
364 
365 namespace {
366 
367 class PhotoCalibSchema {
368 public:
369  table::Schema schema;
370  table::Key<double> calibrationMean;
371  table::Key<double> calibrationErr;
372  table::Key<table::Flag> isConstant;
373  table::Key<int> field;
374  table::Key<int> version;
375 
376  // No copying
377  PhotoCalibSchema(PhotoCalibSchema const &) = delete;
378  PhotoCalibSchema &operator=(PhotoCalibSchema const &) = delete;
379  // No moving
380  PhotoCalibSchema(PhotoCalibSchema &&) = delete;
381  PhotoCalibSchema &operator=(PhotoCalibSchema &&) = delete;
382 
383  static PhotoCalibSchema const &get() {
384  static PhotoCalibSchema const instance;
385  return instance;
386  }
387 
388 private:
389  PhotoCalibSchema()
390  : schema(),
391  calibrationMean(schema.addField<double>(
392  "calibrationMean", "mean calibration on this PhotoCalib's domain", "count")),
394  schema.addField<double>("calibrationErr", "1-sigma error on calibrationMean", "count")),
395  isConstant(schema.addField<table::Flag>("isConstant", "Is this spatially-constant?")),
396  field(schema.addField<int>("field", "archive ID of the BoundedField object")),
397  version(schema.addField<int>("version", "version of this PhotoCalib")) {}
398 };
399 
400 class PhotoCalibFactory : public table::io::PersistableFactory {
401 public:
402  PTR(table::io::Persistable)
403  read(InputArchive const &archive, CatalogVector const &catalogs) const override {
404  table::BaseRecord const &record = catalogs.front().front();
405  PhotoCalibSchema const &keys = PhotoCalibSchema::get();
406  int version = getVersion(record);
407  if (version < 1) {
408  throw(pex::exceptions::RuntimeError("Unsupported version (version 0 was defined in maggies): " +
410  }
411  return std::make_shared<PhotoCalib>(record.get(keys.calibrationMean), record.get(keys.calibrationErr),
412  archive.get<afw::math::BoundedField>(record.get(keys.field)),
413  record.get(keys.isConstant));
414  }
415 
416  PhotoCalibFactory(std::string const &name) : afw::table::io::PersistableFactory(name) {}
417 
418 protected:
419  int getVersion(table::BaseRecord const &record) const {
420  int version = -1;
421  try {
422  std::string versionName = "version";
423  auto versionKey = record.getSchema().find<int>(versionName);
424  version = record.get(versionKey.key);
425  } catch (const pex::exceptions::NotFoundError &) {
426  // un-versioned files are version 0
427  version = 0;
428  }
429  return version;
430  }
431 };
432 
433 std::string getPhotoCalibPersistenceName() { return "PhotoCalib"; }
434 
435 PhotoCalibFactory registration(getPhotoCalibPersistenceName());
436 
437 } // namespace
438 
443 namespace {
444 int const CALIB_TABLE_CURRENT_VERSION = 2; // final version of Calib in ExposureTable
445 std::string const EXPTIME_FIELD_NAME = "exptime"; // name of exposure time field (no longer used)
446 
447 class CalibKeys {
448 public:
449  table::Schema schema;
450  table::Key<std::int64_t> midTime;
451  table::Key<double> expTime;
452  table::Key<double> fluxMag0;
453  table::Key<double> fluxMag0Err;
454 
455  // No copying
456  CalibKeys(const CalibKeys &) = delete;
457  CalibKeys &operator=(const CalibKeys &) = delete;
458 
459  // No moving
460  CalibKeys(CalibKeys &&) = delete;
461  CalibKeys &operator=(CalibKeys &&) = delete;
462 
463  CalibKeys(int tableVersion = CALIB_TABLE_CURRENT_VERSION)
464  : schema(), midTime(), expTime(), fluxMag0(), fluxMag0Err() {
465  if (tableVersion == 1) {
466  // obsolete fields
468  "midtime", "middle of the time of the exposure relative to Unix epoch", "ns");
469  expTime = schema.addField<double>(EXPTIME_FIELD_NAME, "exposure time", "s");
470  }
471  fluxMag0 = schema.addField<double>("fluxmag0", "flux of a zero-magnitude object", "count");
472  fluxMag0Err = schema.addField<double>("fluxmag0.err", "1-sigma error on fluxmag0", "count");
473  }
474 };
475 
476 class CalibFactory : public table::io::PersistableFactory {
477 public:
478  std::shared_ptr<table::io::Persistable> read(InputArchive const &archive,
479  CatalogVector const &catalogs) const override {
480  // table version is not persisted, so we don't have a clean way to determine the version;
481  // the hack is version = 1 if exptime found, else current
482  int tableVersion = 1;
483  try {
484  catalogs.front().getSchema().find<double>(EXPTIME_FIELD_NAME);
485  } catch (pex::exceptions::NotFoundError const&) {
486  tableVersion = CALIB_TABLE_CURRENT_VERSION;
487  }
488 
489  CalibKeys const keys{tableVersion};
490  LSST_ARCHIVE_ASSERT(catalogs.size() == 1u);
491  LSST_ARCHIVE_ASSERT(catalogs.front().size() == 1u);
492  LSST_ARCHIVE_ASSERT(catalogs.front().getSchema() == keys.schema);
493  table::BaseRecord const &record = catalogs.front().front();
494 
495  double calibration = utils::referenceFlux / record.get(keys.fluxMag0);
496  double calibrationErr =
497  utils::referenceFlux * record.get(keys.fluxMag0Err) / std::pow(record.get(keys.fluxMag0), 2);
498  return std::make_shared<PhotoCalib>(calibration, calibrationErr);
499  }
500 
501  explicit CalibFactory(std::string const &name) : table::io::PersistableFactory(name) {}
502 };
503 
504 std::string getCalibPersistenceName() { return "Calib"; }
505 
506 CalibFactory calibRegistration(getCalibPersistenceName());
507 
508 } // namespace
509 
510 std::string PhotoCalib::getPersistenceName() const { return getPhotoCalibPersistenceName(); }
511 
513  PhotoCalibSchema const &keys = PhotoCalibSchema::get();
514  table::BaseCatalog catalog = handle.makeCatalog(keys.schema);
515  auto record = catalog.addNew();
516  record->set(keys.calibrationMean, _calibrationMean);
517  record->set(keys.calibrationErr, _calibrationErr);
518  record->set(keys.isConstant, _isConstant);
519  record->set(keys.field, handle.put(_calibration));
520  record->set(keys.version, SERIALIZATION_VERSION);
521  handle.saveCatalog(catalog);
522 }
523 
524 // ------------------- private/protected helpers -------------------
525 
526 double PhotoCalib::evaluate(lsst::geom::Point<double, 2> const &point) const {
527  if (_isConstant)
528  return _calibrationMean;
529  else
530  return _calibration->evaluate(point);
531 }
532 
533 ndarray::Array<double, 1> PhotoCalib::evaluateArray(ndarray::Array<double, 1> const &xx,
534  ndarray::Array<double, 1> const &yy) const {
535  if (_isConstant) {
536  ndarray::Array<double, 1> result = ndarray::allocate(ndarray::makeVector(xx.size()));
537  result.deep() = _calibrationMean;
538  return result;
539  } else {
540  return _calibration->evaluate(xx, yy);
541  }
542 }
543 
544 ndarray::Array<double, 1> PhotoCalib::evaluateCatalog(afw::table::SourceCatalog const &sourceCatalog) const {
545  ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
546  ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
547  size_t i = 0;
548  for (auto const &rec : sourceCatalog) {
549  auto point = rec.getCentroid();
550  xx[i] = point.getX();
551  yy[i] = point.getY();
552  ++i;
553  }
554  return evaluateArray(xx, yy);
555 }
556 
557 void PhotoCalib::instFluxToNanojanskyArray(afw::table::SourceCatalog const &sourceCatalog,
558  std::string const &instFluxField,
559  ndarray::Array<double, 2, 2> result) const {
560  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
561  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
562 
563  auto calibration = evaluateCatalog(sourceCatalog);
564  int i = 0;
565  auto iter = result.begin();
566  for (auto const &rec : sourceCatalog) {
567  double instFlux = rec.get(instFluxKey);
568  double instFluxErr = rec.get(instFluxErrKey);
569  double nanojansky = toNanojansky(instFlux, calibration[i]);
570  (*iter)[0] = nanojansky;
571  (*iter)[1] = toNanojanskyErr(instFlux, instFluxErr, calibration[i], _calibrationErr, nanojansky);
572  ++iter;
573  ++i;
574  }
575 }
576 
577 void PhotoCalib::instFluxToMagnitudeArray(afw::table::SourceCatalog const &sourceCatalog,
578  std::string const &instFluxField,
579  ndarray::Array<double, 2, 2> result) const {
580  auto instFluxKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFlux").key;
581  auto instFluxErrKey = sourceCatalog.getSchema().find<double>(instFluxField + "_instFluxErr").key;
582 
583  auto calibration = evaluateCatalog(sourceCatalog);
584  auto iter = result.begin();
585  int i = 0;
586  for (auto const &rec : sourceCatalog) {
587  double instFlux = rec.get(instFluxKey);
588  double instFluxErr = rec.get(instFluxErrKey);
589  (*iter)[0] = toMagnitude(instFlux, calibration[i]);
590  (*iter)[1] = toMagnitudeErr(instFlux, instFluxErr, calibration[i], _calibrationErr);
591  ++iter;
592  ++i;
593  }
594 }
595 
597  auto key = "FLUXMAG0";
598  if (metadata.exists(key)) {
599  double instFluxMag0 = metadata.getAsDouble(key);
600  if (strip) metadata.remove(key);
601 
602  double instFluxMag0Err = 0.0;
603  key = "FLUXMAG0ERR";
604  if (metadata.exists(key)) {
605  instFluxMag0Err = metadata.getAsDouble(key);
606  if (strip) metadata.remove(key);
607  }
608  return makePhotoCalibFromCalibZeroPoint(instFluxMag0, instFluxMag0Err);
609  } else {
610  return nullptr;
611  }
612 }
613 
614 std::shared_ptr<PhotoCalib> makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err) {
615  double calibration = utils::referenceFlux / instFluxMag0;
616  double calibrationErr = utils::referenceFlux * instFluxMag0Err / std::pow(instFluxMag0, 2);
617  return std::make_shared<PhotoCalib>(calibration, calibrationErr);
618 }
619 
620 } // namespace image
621 } // namespace afw
622 } // namespace lsst
py::object result
Definition: _schema.cc:430
table::Key< std::string > name
Definition: Amplifier.cc:116
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition: Exception.h:48
Utilities for converting between flux and magnitude in C++.
table::Key< double > calibrationMean
Definition: PhotoCalib.cc:370
table::Schema schema
Definition: PhotoCalib.cc:369
table::Key< std::int64_t > midTime
Definition: PhotoCalib.cc:450
table::Key< double > calibrationErr
Definition: PhotoCalib.cc:371
table::Key< double > fluxMag0
Definition: PhotoCalib.cc:452
table::Key< double > expTime
Definition: PhotoCalib.cc:451
table::Key< table::Flag > isConstant
Definition: PhotoCalib.cc:372
table::Key< double > fluxMag0Err
Definition: PhotoCalib.cc:453
table::Key< int > version
Definition: PhotoCalib.cc:374
table::Key< int > field
Definition: PhotoCalib.cc:373
Implementation of the Photometric Calibration class.
ItemVariant const * other
Definition: Schema.cc:56
Key< U > key
Definition: Schema.cc:281
std::ostream * os
Definition: Schema.cc:746
SchemaMapper * mapper
Definition: SchemaMapper.cc:78
#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
A class to manipulate images, masks, and variance as a single object.
Definition: MaskedImage.h:73
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
Definition: MaskedImage.h:1079
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
Definition: MaskedImage.h:1046
The photometric calibration of an exposure.
Definition: PhotoCalib.h:114
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Definition: BaseRecord.h:80
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition: BaseRecord.h:151
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
Key< T > addField(Field< T > const &field, bool doReplace=false)
Add a new field to the Schema, and return the associated Key.
Definition: Schema.cc:668
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition: Schema.cc:656
A mapping between the keys of two Schemas, used to copy data between them.
Definition: SchemaMapper.h:21
Record class that contains measurements made on a single exposure.
Definition: Source.h:80
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
Definition: Source.h:572
An object passed to Persistable::write to allow it to persist itself.
void saveCatalog(BaseCatalog const &catalog)
Save a catalog in the archive.
BaseCatalog makeCatalog(Schema const &schema)
Return a new, empty catalog with the given schema.
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...
static std::shared_ptr< T > dynamicCast(std::shared_ptr< Persistable > const &ptr)
Dynamically cast a shared_ptr.
Definition: Persistable.cc:18
Interface supporting iteration over heterogenous containers.
Definition: Storable.h:58
Class for storing generic metadata.
Definition: PropertySet.h:67
virtual void remove(std::string const &name)
Remove all values for a property name (possibly hierarchical).
bool exists(std::string const &name) const
Determine if a name (possibly hierarchical) exists.
double getAsDouble(std::string const &name) const
Get the last value for any arithmetic property name (possibly hierarchical).
Reports errors in the logical structure of the program.
Definition: Runtime.h:46
Reports attempts to access elements using an invalid key.
Definition: Runtime.h:151
T emplace_back(T... args)
T find(T... args)
bool strip
Definition: fits.cc:911
T hypot(T... args)
def scale(algorithm, min, max=None, frame=None)
Definition: ds9.py:108
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values.
Definition: PhotoCalib.cc:614
bool operator==(FilterProperty const &rhs) const noexcept
Return true iff two FilterProperties are identical.
std::string getPersistenceName() const override
std::shared_ptr< PhotoCalib > makePhotoCalibFromMetadata(daf::base::PropertySet &metadata, bool strip=false)
Construct a PhotoCalib from FITS FLUXMAG0/FLUXMAG0ERR keywords.
Definition: PhotoCalib.cc:596
void write(OutputArchiveHandle &handle) const override
FilterProperty & operator=(FilterProperty const &)=default
std::ostream & operator<<(std::ostream &os, PhotoCalib const &photoCalib)
Definition: PhotoCalib.cc:257
SortedCatalogT< SourceRecord > SourceCatalog
Definition: fwd.h:85
Angle abs(Angle const &a)
Definition: Angle.h:106
double nanojanskyToABMagnitude(double flux)
Convert a flux in nanojansky to AB magnitude.
Definition: Magnitude.cc:30
double ABMagnitudeToNanojansky(double magnitude)
Convert an AB magnitude to a flux in nanojansky.
Definition: Magnitude.cc:32
const double referenceFlux
The Oke & Gunn (1983) AB magnitude reference flux, in nJy (often approximated as 3631....
Definition: Magnitude.h:46
A base class for image defects.
T pow(T... args)
T setprecision(T... args)
T size(T... args)
T str(T... args)
A value and its error.
Definition: PhotoCalib.h:51
A description of a field in a table.
Definition: Field.h:24
Key< int > photoCalib
Definition: Exposure.cc:67
T to_string(T... args)