42 std::shared_ptr<table::io::Persistable>
const &);
51 s <<
"value=" << measurement.value <<
", error=" << measurement.error;
57int const SERIALIZATION_VERSION = 1;
59double toNanojansky(
double instFlux,
double scale) {
return instFlux * scale; }
61double toMagnitude(
double instFlux,
double scale) {
65double toInstFluxFromMagnitude(
double magnitude,
double scale) {
70double toInstFluxFromNanojansky(
double nanojansky,
double scale) {
return nanojansky /
scale; }
72double toNanojanskyErr(
double instFluxErr,
double scale) {
73 return instFluxErr *
scale;
76double toMagnitudeErr(
double instFlux,
double instFluxErr) {
77 return 2.5 /
std::log(10.0) * (instFluxErr / instFlux);
85 return toNanojansky(instFlux, evaluate(point));
89 return toNanojansky(instFlux, _calibrationMean);
94 double calibration, error, nanojansky;
95 calibration = evaluate(point);
96 nanojansky = toNanojansky(instFlux, calibration);
97 error = toNanojanskyErr(instFluxErr, calibration);
102 double nanojansky = toNanojansky(instFlux, _calibrationMean);
103 double error = toNanojanskyErr(instFluxErr, _calibrationMean);
110 auto instFluxKey = sourceRecord.
getSchema().
find<
double>(instFluxField +
"_instFlux").key;
111 auto instFluxErrKey = sourceRecord.
getSchema().
find<
double>(instFluxField +
"_instFluxErr").key;
116 ndarray::Array<double, 2, 2> result =
117 ndarray::allocate(ndarray::makeVector(
int(sourceCatalog.size()), 2));
118 instFluxToNanojanskyArray(sourceCatalog, instFluxField, result);
124 auto instFluxKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFlux").key;
125 auto instFluxErrKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFluxErr").key;
126 auto nanojanskyKey = sourceCatalog.getSchema().find<
double>(outField +
"_flux").key;
127 auto nanojanskyErrKey = sourceCatalog.getSchema().find<
double>(outField +
"_fluxErr").key;
128 for (
auto &record : sourceCatalog) {
130 record.getCentroid());
131 record.set(nanojanskyKey, result.value);
132 record.set(nanojanskyErrKey, result.error);
139 return toMagnitude(instFlux, evaluate(point));
143 return toMagnitude(instFlux, _calibrationMean);
148 double calibration, error, magnitude;
149 calibration = evaluate(point);
150 magnitude = toMagnitude(instFlux, calibration);
151 error = toMagnitudeErr(instFlux, instFluxErr);
156 double magnitude = toMagnitude(instFlux, _calibrationMean);
157 double error = toMagnitudeErr(instFlux, instFluxErr);
164 auto instFluxKey = sourceRecord.
getSchema().
find<
double>(instFluxField +
"_instFlux").key;
165 auto instFluxErrKey = sourceRecord.
getSchema().
find<
double>(instFluxField +
"_instFluxErr").key;
171 ndarray::Array<double, 2, 2> result =
172 ndarray::allocate(ndarray::makeVector(
int(sourceCatalog.size()), 2));
173 instFluxToMagnitudeArray(sourceCatalog, instFluxField, result);
179 auto instFluxKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFlux").key;
180 auto instFluxErrKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFluxErr").key;
181 auto magKey = sourceCatalog.getSchema().find<
double>(outField +
"_mag").key;
182 auto magErrKey = sourceCatalog.getSchema().find<
double>(outField +
"_magErr").key;
183 for (
auto &record : sourceCatalog) {
185 record.getCentroid());
186 record.set(magKey, result.value);
187 record.set(magErrKey, result.error);
194 return toInstFluxFromMagnitude(magnitude, _calibrationMean);
198 return toInstFluxFromMagnitude(magnitude, evaluate(point));
202 return toInstFluxFromNanojansky(nanojansky, _calibrationMean);
206 return toInstFluxFromNanojansky(nanojansky, evaluate(point));
210 return *(_calibration) / _calibrationMean;
218 return (_calibrationMean == rhs._calibrationMean && _calibrationErr == rhs._calibrationErr &&
219 (*_calibration) == *(rhs._calibration));
223 return calibration->mean();
227 return std::make_unique<PhotoCalib>(*
this);
233 buffer <<
"spatially constant with ";
235 buffer << *_calibration <<
" with ";
236 buffer <<
"mean: " << _calibrationMean <<
" error: " << _calibrationErr;
245 return os << photoCalib.toString();
252 result *= _calibrationMean;
254 _calibration->multiplyImage(result,
true);
261 bool includeScaleUncertainty)
const {
269 result /= _calibrationMean;
271 _calibration->divideImage(result,
true);
277 bool includeScaleUncertainty)
const {
283 auto const &inSchema = catalog.getSchema();
299 keys.reserve(instFluxFields.
size());
300 for (
auto const &field : instFluxFields) {
302 newKey.instFlux = inSchema[inSchema.join(field,
"instFlux")];
304 mapper.
addOutputField(FieldD(inSchema.join(field,
"flux"),
"calibrated flux",
"nJy"),
true);
306 FieldD(inSchema.join(field,
"mag"),
"calibrated magnitude",
"mag(AB)"),
true);
308 newKey.instFluxErr = inSchema.find<
double>(inSchema.join(field,
"instFluxErr")).key;
310 FieldD(inSchema.join(field,
"fluxErr"),
"calibrated flux uncertainty",
"nJy"),
true);
312 FieldD(inSchema.join(field,
"magErr"),
"calibrated magnitude uncertainty",
"mag(AB)"),
317 keys.emplace_back(newKey);
322 output.insert(mapper, output.begin(), catalog.begin(), catalog.end());
324 auto calibration = evaluateCatalog(output);
328 for (
auto &rec : output) {
329 for (
auto &key : keys) {
330 double instFlux = rec.get(key.instFlux);
331 double nanojansky = toNanojansky(instFlux, calibration[iRec]);
332 rec.set(key.flux, nanojansky);
333 rec.set(key.mag, toMagnitude(instFlux, calibration[iRec]));
334 if (key.instFluxErr.isValid()) {
335 double instFluxErr = rec.get(key.instFluxErr);
336 rec.set(key.fluxErr, toNanojanskyErr(instFluxErr, calibration[iRec]));
337 rec.set(key.magErr, toMagnitudeErr(instFlux, instFluxErr));
349 for (
auto const &name : catalog.getSchema().getNames()) {
351 if (name.size() > SUFFIX.
size() + 1 &&
352 name.compare(name.size() - SUFFIX.
size(), SUFFIX.
size(), SUFFIX) == 0) {
353 instFluxFields.
emplace_back(name.substr(0, name.size() - 9));
363class PhotoCalibSchema {
373 PhotoCalibSchema(PhotoCalibSchema
const &) =
delete;
374 PhotoCalibSchema &operator=(PhotoCalibSchema
const &) =
delete;
376 PhotoCalibSchema(PhotoCalibSchema &&) =
delete;
377 PhotoCalibSchema &operator=(PhotoCalibSchema &&) =
delete;
379 static PhotoCalibSchema
const &get() {
380 static PhotoCalibSchema
const instance;
387 calibrationMean(schema.addField<double>(
388 "calibrationMean",
"mean calibration on this PhotoCalib's domain",
"count")),
390 schema.addField<double>(
"calibrationErr",
"1-sigma error on calibrationMean",
"count")),
391 isConstant(schema.addField<
table::Flag>(
"isConstant",
"Is this spatially-constant?")),
392 field(schema.addField<int>(
"field",
"archive ID of the BoundedField object")),
393 version(schema.addField<int>(
"version",
"version of this PhotoCalib")) {}
396class PhotoCalibFactory :
public table::io::PersistableFactory {
399 CatalogVector
const &catalogs)
const override {
400 table::BaseRecord
const &record = catalogs.front().front();
401 PhotoCalibSchema
const &keys = PhotoCalibSchema::get();
402 int version = getVersion(record);
404 throw(pex::exceptions::RuntimeError(
"Unsupported version (version 0 was defined in maggies): " +
408 archive.get<afw::math::BoundedField>(record.get(
keys.field)),
409 record.get(
keys.isConstant));
412 PhotoCalibFactory(std::string
const &name) : afw::
table::io::PersistableFactory(
name) {}
418 std::string versionName =
"version";
419 auto versionKey = record.getSchema().find<
int>(versionName);
420 version = record.get(versionKey.key);
421 }
catch (
const pex::exceptions::NotFoundError &) {
429std::string getPhotoCalibPersistenceName() {
return "PhotoCalib"; }
431PhotoCalibFactory registration(getPhotoCalibPersistenceName());
440int const CALIB_TABLE_CURRENT_VERSION = 2;
441std::string
const EXPTIME_FIELD_NAME =
"exptime";
452 CalibKeys(
const CalibKeys &) =
delete;
453 CalibKeys &operator=(
const CalibKeys &) =
delete;
456 CalibKeys(CalibKeys &&) =
delete;
457 CalibKeys &operator=(CalibKeys &&) =
delete;
459 CalibKeys(
int tableVersion = CALIB_TABLE_CURRENT_VERSION)
460 :
schema(), midTime(), expTime(), fluxMag0(), fluxMag0Err() {
461 if (tableVersion == 1) {
463 midTime =
schema.addField<std::int64_t>(
464 "midtime",
"middle of the time of the exposure relative to Unix epoch",
"ns");
465 expTime =
schema.addField<
double>(EXPTIME_FIELD_NAME,
"exposure time",
"s");
467 fluxMag0 =
schema.addField<
double>(
"fluxmag0",
"flux of a zero-magnitude object",
"count");
468 fluxMag0Err =
schema.addField<
double>(
"fluxmag0.err",
"1-sigma error on fluxmag0",
"count");
474 std::shared_ptr<table::io::Persistable>
read(InputArchive
const &archive,
475 CatalogVector
const &catalogs)
const override {
478 int tableVersion = 1;
480 catalogs.front().getSchema().find<
double>(EXPTIME_FIELD_NAME);
481 }
catch (pex::exceptions::NotFoundError
const &) {
482 tableVersion = CALIB_TABLE_CURRENT_VERSION;
485 CalibKeys
const keys{tableVersion};
497 explicit CalibFactory(std::string
const &name) :
table::io::PersistableFactory(
name) {}
500std::string getCalibPersistenceName() {
return "Calib"; }
502CalibFactory calibRegistration(getCalibPersistenceName());
509 PhotoCalibSchema
const &keys = PhotoCalibSchema::get();
511 auto record = catalog.
addNew();
512 record->set(keys.calibrationMean, _calibrationMean);
513 record->set(keys.calibrationErr, _calibrationErr);
514 record->set(keys.isConstant, _isConstant);
515 record->set(keys.field, handle.
put(_calibration));
516 record->set(keys.version, SERIALIZATION_VERSION);
524 return _calibrationMean;
526 return _calibration->evaluate(point);
529ndarray::Array<double, 1> PhotoCalib::evaluateArray(ndarray::Array<double, 1>
const &xx,
530 ndarray::Array<double, 1>
const &yy)
const {
532 ndarray::Array<double, 1> result = ndarray::allocate(ndarray::makeVector(xx.size()));
533 result.deep() = _calibrationMean;
536 return _calibration->evaluate(xx, yy);
541 ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
542 ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
544 for (
auto const &rec : sourceCatalog) {
545 auto point = rec.getCentroid();
546 xx[i] = point.getX();
547 yy[i] = point.getY();
550 return evaluateArray(xx, yy);
554 std::string
const &instFluxField,
555 ndarray::Array<double, 2, 2> result)
const {
556 auto instFluxKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFlux").key;
557 auto instFluxErrKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFluxErr").key;
559 auto calibration = evaluateCatalog(sourceCatalog);
562 for (
auto const &rec : sourceCatalog) {
563 double instFlux = rec.get(instFluxKey);
564 double instFluxErr = rec.get(instFluxErrKey);
565 double nanojansky = toNanojansky(instFlux, calibration[i]);
566 (*iter)[0] = nanojansky;
567 (*iter)[1] = toNanojanskyErr(instFluxErr, calibration[i]);
574 std::string
const &instFluxField,
575 ndarray::Array<double, 2, 2> result)
const {
576 auto instFluxKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFlux").key;
577 auto instFluxErrKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFluxErr").key;
579 auto calibration = evaluateCatalog(sourceCatalog);
582 for (
auto const &rec : sourceCatalog) {
583 double instFlux = rec.get(instFluxKey);
584 double instFluxErr = rec.get(instFluxErrKey);
585 (*iter)[0] = toMagnitude(instFlux, calibration[i]);
586 (*iter)[1] = toMagnitudeErr(instFlux, instFluxErr);
593 auto key =
"FLUXMAG0";
594 if (metadata.
exists(key)) {
596 if (strip) metadata.
remove(key);
598 double instFluxMag0Err = 0.0;
600 if (metadata.
exists(key)) {
602 if (strip) metadata.
remove(key);
#define LSST_EXCEPT(type,...)
Create an exception with a given type.
Utilities for converting between flux and magnitude in C++.
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.
afw::table::SourceCatalog calibrateCatalog(afw::table::SourceCatalog const &catalog, std::vector< std::string > const &instFluxFields) const
Return a flux calibrated catalog, with new _flux, _fluxErr, _mag, and _magErr fields.
bool operator==(PhotoCalib const &rhs) const
Two PhotoCalibs are equal if their component bounded fields and calibrationErr are equal.
double instFluxToNanojansky(double instFlux, lsst::geom::Point< double, 2 > const &point) const
Convert instFlux in ADU to nJy at a point in the BoundedField.
std::string getPersistenceName() const override
Return the unique name used to persist this object and look up its factory.
PhotoCalib(PhotoCalib const &)=default
double instFluxToMagnitude(double instFlux, lsst::geom::Point< double, 2 > const &point) const
Convert instFlux in ADU to AB magnitude.
std::string toString() const override
Create a string representation of this object.
bool equals(typehandling::Storable const &other) const noexcept override
Compare this object to another Storable.
MaskedImage< float > uncalibrateImage(MaskedImage< float > const &maskedImage) const
Return a un-calibrated image, with pixel values in ADU (or whatever the original input to this photoC...
std::shared_ptr< typehandling::Storable > cloneStorable() const override
Create a new PhotoCalib that is a copy of this one.
double nanojanskyToInstFlux(double nanojansky, lsst::geom::Point< double, 2 > const &point) const
Convert nanojansky to instFlux (ADU).
std::shared_ptr< afw::math::BoundedField > computeScalingTo(std::shared_ptr< PhotoCalib > other) const
Calculates the scaling between this PhotoCalib and another PhotoCalib.
void write(OutputArchiveHandle &handle) const override
Write the object to one or more catalogs.
std::shared_ptr< afw::math::BoundedField > computeScaledCalibration() const
Calculates the spatially-variable calibration, normalized by the mean in the valid domain.
MaskedImage< float > calibrateImage(MaskedImage< float > const &maskedImage) const
Return a flux calibrated image, with pixel values in nJy.
double magnitudeToInstFlux(double magnitude, lsst::geom::Point< double, 2 > const &point) const
Convert AB magnitude to instFlux (ADU).
Base class for all records.
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.
Defines the fields and offsets for 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.
Schema const getOutputSchema() const
Return the output schema (copy-on-write).
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.
void addMinimalSchema(Schema const &minimal, bool doMap=true)
Add the given minimal schema to the output schema.
Record class that contains measurements made on a single exposure.
CentroidSlotDefinition::MeasValue getCentroid() const
Get the value of the Centroid slot measurement.
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.
A base class for factory classes used to reconstruct objects from records.
io::OutputArchiveHandle OutputArchiveHandle
Interface supporting iteration over heterogenous containers.
static bool singleClassEquals(T const &lhs, Storable const &rhs)
Test if a Storable is of a particular class and equal to another object.
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.
T emplace_back(T... args)
scale(algorithm, min, max=None, frame=None)
std::shared_ptr< PhotoCalib > makePhotoCalibFromCalibZeroPoint(double instFluxMag0, double instFluxMag0Err)
Construct a PhotoCalib from the deprecated Calib-style instFluxMag0/instFluxMag0Err values.
std::shared_ptr< PhotoCalib > makePhotoCalibFromMetadata(daf::base::PropertySet &metadata, bool strip=false)
Construct a PhotoCalib from FITS FLUXMAG0/FLUXMAG0ERR keywords.
std::ostream & operator<<(std::ostream &os, Measurement const &measurement)
CatalogT< BaseRecord > BaseCatalog
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