51 s <<
"value=" << measurement.value <<
", error=" << measurement.error;
57int const SERIALIZATION_VERSION = 1;
59double toNanojansky(
double instFlux,
double scale) {
return instFlux * scale; }
63double toInstFluxFromMagnitude(
double magnitude,
double scale) {
68double toNanojanskyErr(
double instFlux,
double instFluxErr,
double scale,
double scaleErr,
70 return std::abs(nanojansky) *
hypot(instFluxErr / instFlux, scaleErr / scale);
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());
97double toMagnitudeErr(
double instFlux,
double instFluxErr,
double scale,
double scaleErr) {
98 return 2.5 /
std::log(10.0) *
hypot(instFluxErr / instFlux, scaleErr / scale);
106 return toNanojansky(instFlux, evaluate(point));
109double PhotoCalib::instFluxToNanojansky(
double instFlux)
const {
110 return toNanojansky(instFlux, _calibrationMean);
109double PhotoCalib::instFluxToNanojansky(
double instFlux)
const {
…}
113Measurement PhotoCalib::instFluxToNanojansky(
double instFlux,
double instFluxErr,
115 double calibration, error, nanojansky;
116 calibration = evaluate(point);
117 nanojansky = toNanojansky(instFlux, calibration);
118 error = toNanojanskyErr(instFlux, instFluxErr, calibration, _calibrationErr, nanojansky);
113Measurement PhotoCalib::instFluxToNanojansky(
double instFlux,
double instFluxErr, {
…}
122Measurement PhotoCalib::instFluxToNanojansky(
double instFlux,
double instFluxErr)
const {
123 double nanojansky = toNanojansky(instFlux, _calibrationMean);
124 double error = toNanojanskyErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr, nanojansky);
122Measurement PhotoCalib::instFluxToNanojansky(
double instFlux,
double instFluxErr)
const {
…}
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);
137 ndarray::Array<double, 2, 2>
result =
138 ndarray::allocate(ndarray::makeVector(
int(sourceCatalog.size()), 2));
139 instFluxToNanojanskyArray(sourceCatalog, instFluxField,
result);
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);
160 return toMagnitude(instFlux, evaluate(point));
163double PhotoCalib::instFluxToMagnitude(
double instFlux)
const {
164 return toMagnitude(instFlux, _calibrationMean);
163double PhotoCalib::instFluxToMagnitude(
double instFlux)
const {
…}
167Measurement PhotoCalib::instFluxToMagnitude(
double instFlux,
double instFluxErr,
169 double calibration, error, magnitude;
170 calibration = evaluate(point);
171 magnitude = toMagnitude(instFlux, calibration);
172 error = toMagnitudeErr(instFlux, instFluxErr, calibration, _calibrationErr);
167Measurement PhotoCalib::instFluxToMagnitude(
double instFlux,
double instFluxErr, {
…}
176Measurement PhotoCalib::instFluxToMagnitude(
double instFlux,
double instFluxErr)
const {
177 double magnitude = toMagnitude(instFlux, _calibrationMean);
178 double error = toMagnitudeErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr);
176Measurement PhotoCalib::instFluxToMagnitude(
double instFlux,
double instFluxErr)
const {
…}
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);
192 ndarray::Array<double, 2, 2>
result =
193 ndarray::allocate(ndarray::makeVector(
int(sourceCatalog.size()), 2));
194 instFluxToMagnitudeArray(sourceCatalog, instFluxField,
result);
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);
214double PhotoCalib::magnitudeToInstFlux(
double magnitude)
const {
215 return toInstFluxFromMagnitude(magnitude, _calibrationMean);
214double PhotoCalib::magnitudeToInstFlux(
double magnitude)
const {
…}
219 return toInstFluxFromMagnitude(magnitude, evaluate(point));
223 return *(_calibration) / _calibrationMean;
231 return (_calibrationMean == rhs._calibrationMean && _calibrationErr == rhs._calibrationErr &&
232 (*_calibration) == *(rhs._calibration));
236 return calibration->mean();
240 return std::make_unique<PhotoCalib>(*
this);
246 buffer <<
"spatially constant with ";
248 buffer << *_calibration <<
" with ";
249 buffer <<
"mean: " << _calibrationMean <<
" error: " << _calibrationErr;
254 return singleClassEquals(*
this, other);
258 return os << photoCalib.toString();
262 bool includeScaleUncertainty)
const {
267 *(
result.getImage()) *= _calibrationMean;
269 _calibration->multiplyImage(*(
result.getImage()),
true);
271 if (includeScaleUncertainty) {
272 toNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(),
273 _calibrationErr,
result.getImage()->getArray(),
274 result.getVariance()->getArray());
276 toNanojanskyVariance(maskedImage.getImage()->getArray(), maskedImage.getVariance()->getArray(), 0,
277 result.getImage()->getArray(),
result.getVariance()->getArray());
285 auto const &inSchema = catalog.getSchema();
287 mapper.addMinimalSchema(inSchema);
301 keys.reserve(instFluxFields.
size());
302 for (
auto const &field : instFluxFields) {
304 newKey.instFlux = inSchema[inSchema.join(field,
"instFlux")];
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);
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)"),
319 keys.emplace_back(newKey);
324 output.insert(
mapper, output.begin(), catalog.begin(), catalog.end());
326 auto calibration = evaluateCatalog(output);
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));
341 toMagnitudeErr(instFlux, instFluxErr, calibration[iRec], _calibrationErr));
353 for (
auto const &name : catalog.getSchema().getNames()) {
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));
360 return calibrateCatalog(catalog, instFluxFields);
367class PhotoCalibSchema {
377 PhotoCalibSchema(PhotoCalibSchema
const &) =
delete;
378 PhotoCalibSchema &operator=(PhotoCalibSchema
const &) =
delete;
380 PhotoCalibSchema(PhotoCalibSchema &&) =
delete;
381 PhotoCalibSchema &operator=(PhotoCalibSchema &&) =
delete;
383 static PhotoCalibSchema
const &get() {
384 static PhotoCalibSchema
const instance;
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")) {}
400class PhotoCalibFactory :
public table::io::PersistableFactory {
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);
408 throw(pex::exceptions::RuntimeError(
"Unsupported version (version 0 was defined in maggies): " +
412 archive.get<afw::math::BoundedField>(record.get(
keys.field)),
413 record.get(
keys.isConstant));
416 PhotoCalibFactory(
std::string const &name) :
afw::table::io::PersistableFactory(name) {}
419 int getVersion(table::BaseRecord
const &record)
const {
423 auto versionKey = record.getSchema().
find<
int>(versionName);
424 version = record.get(versionKey.key);
425 }
catch (
const pex::exceptions::NotFoundError &) {
433std::string getPhotoCalibPersistenceName() {
return "PhotoCalib"; }
435PhotoCalibFactory registration(getPhotoCalibPersistenceName());
444int const CALIB_TABLE_CURRENT_VERSION = 2;
456 CalibKeys(
const CalibKeys &) =
delete;
457 CalibKeys &operator=(
const CalibKeys &) =
delete;
460 CalibKeys(CalibKeys &&) =
delete;
461 CalibKeys &operator=(CalibKeys &&) =
delete;
463 CalibKeys(
int tableVersion = CALIB_TABLE_CURRENT_VERSION)
465 if (tableVersion == 1) {
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");
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");
476class CalibFactory :
public table::io::PersistableFactory {
479 CatalogVector
const &catalogs)
const override {
482 int tableVersion = 1;
484 catalogs.front().getSchema().find<
double>(EXPTIME_FIELD_NAME);
485 }
catch (pex::exceptions::NotFoundError
const&) {
486 tableVersion = CALIB_TABLE_CURRENT_VERSION;
489 CalibKeys
const keys{tableVersion};
493 table::BaseRecord
const &record = catalogs.front().front();
501 explicit CalibFactory(
std::string const &name) : table::io::PersistableFactory(name) {}
504std::string getCalibPersistenceName() {
return "Calib"; }
506CalibFactory calibRegistration(getCalibPersistenceName());
510std::string PhotoCalib::getPersistenceName()
const {
return getPhotoCalibPersistenceName(); }
513 PhotoCalibSchema
const &keys = PhotoCalibSchema::get();
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);
528 return _calibrationMean;
530 return _calibration->evaluate(point);
533ndarray::Array<double, 1> PhotoCalib::evaluateArray(ndarray::Array<double, 1>
const &xx,
534 ndarray::Array<double, 1>
const &yy)
const {
536 ndarray::Array<double, 1>
result = ndarray::allocate(ndarray::makeVector(xx.size()));
537 result.deep() = _calibrationMean;
540 return _calibration->evaluate(xx, yy);
545 ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
546 ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
548 for (
auto const &rec : sourceCatalog) {
549 auto point = rec.getCentroid();
550 xx[i] = point.getX();
551 yy[i] = point.getY();
554 return evaluateArray(xx, yy);
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;
563 auto calibration = evaluateCatalog(sourceCatalog);
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);
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;
583 auto calibration = evaluateCatalog(sourceCatalog);
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);
597 auto key =
"FLUXMAG0";
598 if (metadata.
exists(key)) {
602 double instFluxMag0Err = 0.0;
604 if (metadata.
exists(key)) {
608 return makePhotoCalibFromCalibZeroPoint(instFluxMag0, instFluxMag0Err);
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
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.
#define LSST_ARCHIVE_ASSERT(EXPR)
An assertion macro used to validate the structure of an InputArchive.
A class to manipulate images, masks, and variance as a single object.
The photometric calibration of an exposure.
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
A class used as a handle to a particular field in a table.
SchemaItem< T > find(std::string const &name) const
Find a SchemaItem in the Schema by name.
A mapping between the keys of two Schemas, used to copy data between them.
Record class that contains measurements made on a single exposure.
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
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.
Class for storing generic metadata.
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.
Reports errors in the logical structure of the program.
Reports attempts to access elements using an invalid key.
Utilities for converting between flux and magnitude in C++.
T emplace_back(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
const double referenceFlux
The Oke & Gunn (1983) AB magnitude reference flux, in nJy (often approximated as 3631....
double ABMagnitudeToNanojansky(double magnitude)
Convert an AB magnitude to a flux in nanojansky.
double nanojanskyToABMagnitude(double flux)
Convert a flux in nanojansky to AB magnitude.
T setprecision(T... args)
A description of a field in a table.
std::shared_ptr< table::io::Persistable > read(table::io::InputArchive const &archive, table::io::CatalogVector const &catalogs) const override