LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
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#include <iomanip>
26
27#include "lsst/geom/Point.h"
33#include "lsst/pex/exceptions.h"
34#include "ndarray.h"
37
38namespace lsst {
39namespace afw {
40
43
44namespace image {
45
46// ------------------- helpers -------------------
47
48std::ostream &operator<<(std::ostream &os, Measurement const &measurement) {
51 s << "value=" << measurement.value << ", error=" << measurement.error;
52 return os << s.str();
53}
54
55namespace {
56
57int const SERIALIZATION_VERSION = 1;
58
59double toNanojansky(double instFlux, double scale) { return instFlux * scale; }
60
61double toMagnitude(double instFlux, double scale) { return utils::nanojanskyToABMagnitude(instFlux * scale); }
62
63double toInstFluxFromMagnitude(double magnitude, double scale) {
64 // Note: flux[nJy] / scale = instFlux[counts]
65 return utils::ABMagnitudeToNanojansky(magnitude) / scale;
66}
67
68double toNanojanskyErr(double instFlux, double instFluxErr, double scale, double scaleErr,
69 double nanojansky) {
70 return std::abs(nanojansky) * hypot(instFluxErr / instFlux, scaleErr / scale);
71}
72
86void 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
97double toMagnitudeErr(double instFlux, double instFluxErr, double scale, double scaleErr) {
98 return 2.5 / std::log(10.0) * hypot(instFluxErr / instFlux, scaleErr / scale);
99}
100
101} // anonymous namespace
102
103// ------------------- Conversions to nanojansky -------------------
104
105double PhotoCalib::instFluxToNanojansky(double instFlux, lsst::geom::Point<double, 2> const &point) const {
106 return toNanojansky(instFlux, evaluate(point));
107}
108
109double PhotoCalib::instFluxToNanojansky(double instFlux) const {
110 return toNanojansky(instFlux, _calibrationMean);
111}
112
113Measurement 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
122Measurement 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
128Measurement 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}
135ndarray::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
143void 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
159double PhotoCalib::instFluxToMagnitude(double instFlux, lsst::geom::Point<double, 2> const &point) const {
160 return toMagnitude(instFlux, evaluate(point));
161}
162
163double PhotoCalib::instFluxToMagnitude(double instFlux) const {
164 return toMagnitude(instFlux, _calibrationMean);
165}
166
167Measurement 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
176Measurement 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
182Measurement 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
190ndarray::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
198void 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
214double PhotoCalib::magnitudeToInstFlux(double magnitude) const {
215 return toInstFluxFromMagnitude(magnitude, _calibrationMean);
216}
217
218double PhotoCalib::magnitudeToInstFlux(double magnitude, lsst::geom::Point<double, 2> const &point) const {
219 return toInstFluxFromMagnitude(magnitude, evaluate(point));
220}
221
222std::shared_ptr<math::BoundedField> PhotoCalib::computeScaledCalibration() const {
223 return *(_calibration) / _calibrationMean;
224}
225
226std::shared_ptr<math::BoundedField> PhotoCalib::computeScalingTo(std::shared_ptr<PhotoCalib> other) const {
227 throw LSST_EXCEPT(pex::exceptions::LogicError, "Not Implemented: See DM-10154.");
228}
229
230bool PhotoCalib::operator==(PhotoCalib const &rhs) const {
231 return (_calibrationMean == rhs._calibrationMean && _calibrationErr == rhs._calibrationErr &&
232 (*_calibration) == *(rhs._calibration));
233}
234
235double PhotoCalib::computeCalibrationMean(std::shared_ptr<afw::math::BoundedField> calibration) const {
236 return calibration->mean();
237}
238
239std::shared_ptr<typehandling::Storable> PhotoCalib::cloneStorable() const {
240 return std::make_unique<PhotoCalib>(*this);
241}
242
243std::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
253bool PhotoCalib::equals(typehandling::Storable const &other) const noexcept {
254 return singleClassEquals(*this, other);
255}
256
257std::ostream &operator<<(std::ostream &os, PhotoCalib const &photoCalib) {
258 return os << photoCalib.toString();
259}
260
261MaskedImage<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
283afw::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;
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);
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
350afw::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
365namespace {
366
367class PhotoCalibSchema {
368public:
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
388private:
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
400class PhotoCalibFactory : public table::io::PersistableFactory {
401public:
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
418protected:
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
433std::string getPhotoCalibPersistenceName() { return "PhotoCalib"; }
434
435PhotoCalibFactory registration(getPhotoCalibPersistenceName());
436
437} // namespace
438
439/*
440 * Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
441 */
442
443namespace {
444int const CALIB_TABLE_CURRENT_VERSION = 2; // final version of Calib in ExposureTable
445std::string const EXPTIME_FIELD_NAME = "exptime"; // name of exposure time field (no longer used)
446
447class CalibKeys {
448public:
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
467 midTime = schema.addField<std::int64_t>(
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
476class CalibFactory : public table::io::PersistableFactory {
477public:
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
504std::string getCalibPersistenceName() { return "Calib"; }
505
506CalibFactory calibRegistration(getCalibPersistenceName());
507
508} // namespace
509
510std::string PhotoCalib::getPersistenceName() const { return getPhotoCalibPersistenceName(); }
511
512void PhotoCalib::write(OutputArchiveHandle &handle) const {
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
526double 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
533ndarray::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
544ndarray::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
557void 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
577void 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
614std::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:429
table::Key< int > field
Definition ApCorrMap.cc:77
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Definition Exception.h:48
table::Key< double > calibrationMean
table::Key< std::int64_t > midTime
table::Key< double > calibrationErr
table::Key< double > fluxMag0
table::Key< double > expTime
table::Key< table::Flag > isConstant
table::Key< double > fluxMag0Err
Implementation of the Photometric Calibration class.
SchemaMapper * mapper
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
Definition Persistable.h:48
table::Schema schema
Definition python.h:134
A class to manipulate images, masks, and variance as a single object.
Definition MaskedImage.h:74
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
Tag types used to declare specialized field types.
Definition misc.h:31
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Definition BaseRecord.h:151
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Definition BaseRecord.h:80
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:489
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
Definition Schema.cc:467
A mapping between the keys of two Schemas, used to copy data between them.
Record class that contains measurements made on a single exposure.
Definition Source.h:78
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
Definition Source.h:570
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.
Interface supporting iteration over heterogenous containers.
Definition Storable.h:58
Class for storing generic metadata.
Definition PropertySet.h:66
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).
A coordinate class intended to represent absolute positions.
Definition Point.h:169
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:930
T hypot(T... args)
T log(T... args)
std::shared_ptr< PhotoCalib > makePhotoCalibFromMetadata(daf::base::PropertySet &metadata, bool strip=false)
Construct a PhotoCalib from FITS FLUXMAG0/FLUXMAG0ERR keywords.
SortedCatalogT< SourceRecord > SourceCatalog
Definition fwd.h:85
const double referenceFlux
The Oke & Gunn (1983) AB magnitude reference flux, in nJy (often approximated as 3631....
Definition Magnitude.h:46
double ABMagnitudeToNanojansky(double magnitude)
Convert an AB magnitude to a flux in nanojansky.
Definition Magnitude.cc:32
double nanojanskyToABMagnitude(double flux)
Convert a flux in nanojansky to AB magnitude.
Definition Magnitude.cc:30
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
T to_string(T... args)
Utilities for converting between flux and magnitude in C++.
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override