50 s <<
"value=" << measurement.
value <<
", error=" << measurement.
error;
56 int const SERIALIZATION_VERSION = 1;
58 double toNanojansky(
double instFlux,
double scale) {
return instFlux *
scale; }
62 double toInstFluxFromMagnitude(
double magnitude,
double scale) {
67 double toNanojanskyErr(
double instFlux,
double instFluxErr,
double scale,
double scaleErr,
85 void toNanojanskyVariance(ndarray::Array<float const, 2, 1>
const &instFlux,
86 ndarray::Array<float const, 2, 1>
const &instFluxVar,
float scaleErr,
87 ndarray::Array<float const, 2, 1>
const &flux, ndarray::Array<float, 2, 1> out) {
88 auto eigenFlux = ndarray::asEigen<Eigen::ArrayXpr>(flux);
89 auto eigenInstFluxVar = ndarray::asEigen<Eigen::ArrayXpr>(instFluxVar);
90 auto eigenInstFlux = ndarray::asEigen<Eigen::ArrayXpr>(instFlux);
91 auto eigenOut = ndarray::asEigen<Eigen::ArrayXpr>(out);
92 eigenOut = eigenFlux.square() *
93 (eigenInstFluxVar / eigenInstFlux.square() + (scaleErr / eigenFlux * eigenInstFlux).square());
96 double toMagnitudeErr(
double instFlux,
double instFluxErr,
double scale,
double scaleErr) {
105 return toNanojansky(instFlux, evaluate(point));
108 double PhotoCalib::instFluxToNanojansky(
double instFlux)
const {
109 return toNanojansky(instFlux, _calibrationMean);
112 Measurement PhotoCalib::instFluxToNanojansky(
double instFlux,
double instFluxErr,
114 double calibration,
error, nanojansky;
115 calibration = evaluate(point);
116 nanojansky = toNanojansky(instFlux, calibration);
117 error = toNanojanskyErr(instFlux, instFluxErr, calibration, _calibrationErr, nanojansky);
121 Measurement PhotoCalib::instFluxToNanojansky(
double instFlux,
double instFluxErr)
const {
122 double nanojansky = toNanojansky(instFlux, _calibrationMean);
123 double error = toNanojanskyErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr, nanojansky);
130 auto instFluxKey = sourceRecord.
getSchema().
find<
double>(instFluxField +
"_instFlux").key;
131 auto instFluxErrKey = sourceRecord.
getSchema().
find<
double>(instFluxField +
"_instFluxErr").key;
132 return instFluxToNanojansky(sourceRecord.
get(instFluxKey), sourceRecord.
get(instFluxErrKey), position);
136 ndarray::Array<double, 2, 2>
result =
137 ndarray::allocate(ndarray::makeVector(
int(sourceCatalog.size()), 2));
138 instFluxToNanojanskyArray(sourceCatalog, instFluxField,
result);
144 auto instFluxKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFlux").key;
145 auto instFluxErrKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFluxErr").key;
146 auto nanojanskyKey = sourceCatalog.getSchema().find<
double>(outField +
"_flux").key;
147 auto nanojanskyErrKey = sourceCatalog.getSchema().find<
double>(outField +
"_fluxErr").key;
148 for (
auto &record : sourceCatalog) {
149 auto result = instFluxToNanojansky(record.get(instFluxKey), record.get(instFluxErrKey),
150 record.getCentroid());
151 record.set(nanojanskyKey,
result.value);
152 record.set(nanojanskyErrKey,
result.error);
159 return toMagnitude(instFlux, evaluate(point));
162 double PhotoCalib::instFluxToMagnitude(
double instFlux)
const {
163 return toMagnitude(instFlux, _calibrationMean);
166 Measurement PhotoCalib::instFluxToMagnitude(
double instFlux,
double instFluxErr,
168 double calibration,
error, magnitude;
169 calibration = evaluate(point);
170 magnitude = toMagnitude(instFlux, calibration);
171 error = toMagnitudeErr(instFlux, instFluxErr, calibration, _calibrationErr);
175 Measurement PhotoCalib::instFluxToMagnitude(
double instFlux,
double instFluxErr)
const {
176 double magnitude = toMagnitude(instFlux, _calibrationMean);
177 double error = toMagnitudeErr(instFlux, instFluxErr, _calibrationMean, _calibrationErr);
184 auto instFluxKey = sourceRecord.
getSchema().
find<
double>(instFluxField +
"_instFlux").key;
185 auto instFluxErrKey = sourceRecord.
getSchema().
find<
double>(instFluxField +
"_instFluxErr").key;
186 return instFluxToMagnitude(sourceRecord.
get(instFluxKey), sourceRecord.
get(instFluxErrKey), position);
191 ndarray::Array<double, 2, 2>
result =
192 ndarray::allocate(ndarray::makeVector(
int(sourceCatalog.size()), 2));
193 instFluxToMagnitudeArray(sourceCatalog, instFluxField,
result);
199 auto instFluxKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFlux").key;
200 auto instFluxErrKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFluxErr").key;
201 auto magKey = sourceCatalog.getSchema().find<
double>(outField +
"_mag").key;
202 auto magErrKey = sourceCatalog.getSchema().find<
double>(outField +
"_magErr").key;
203 for (
auto &record : sourceCatalog) {
204 auto result = instFluxToMagnitude(record.get(instFluxKey), record.get(instFluxErrKey),
205 record.getCentroid());
206 record.set(magKey,
result.value);
207 record.set(magErrKey,
result.error);
213 double PhotoCalib::magnitudeToInstFlux(
double magnitude)
const {
214 return toInstFluxFromMagnitude(magnitude, _calibrationMean);
218 return toInstFluxFromMagnitude(magnitude, evaluate(point));
222 return *(_calibration) / _calibrationMean;
230 return (_calibrationMean == rhs._calibrationMean && _calibrationErr == rhs._calibrationErr &&
231 (*_calibration) == *(rhs._calibration));
235 return calibration->mean();
239 return std::make_unique<PhotoCalib>(*
this);
245 buffer <<
"spatially constant with ";
247 buffer << *_calibration <<
" with ";
248 buffer <<
"mean: " << _calibrationMean <<
" error: " << _calibrationErr;
253 return singleClassEquals(*
this, other);
261 bool includeScaleUncertainty)
const {
266 *(
result.getImage()) *= _calibrationMean;
268 _calibration->multiplyImage(*(
result.getImage()),
true);
270 if (includeScaleUncertainty) {
271 toNanojanskyVariance(maskedImage.
getImage()->getArray(), maskedImage.
getVariance()->getArray(),
272 _calibrationErr,
result.getImage()->getArray(),
273 result.getVariance()->getArray());
275 toNanojanskyVariance(maskedImage.
getImage()->getArray(), maskedImage.
getVariance()->getArray(), 0,
276 result.getImage()->getArray(),
result.getVariance()->getArray());
284 auto const &inSchema = catalog.getSchema();
286 mapper.addMinimalSchema(inSchema);
300 keys.reserve(instFluxFields.
size());
301 for (
auto const &
field : instFluxFields) {
303 newKey.instFlux = inSchema[inSchema.join(
field,
"instFlux")];
305 mapper.addOutputField(FieldD(inSchema.join(
field,
"flux"),
"calibrated flux",
"nJy"),
true);
306 newKey.mag =
mapper.addOutputField(
307 FieldD(inSchema.join(
field,
"mag"),
"calibrated magnitude",
"mag(AB)"),
true);
309 newKey.instFluxErr = inSchema.find<
double>(inSchema.join(
field,
"instFluxErr")).key;
310 newKey.fluxErr =
mapper.addOutputField(
311 FieldD(inSchema.join(
field,
"fluxErr"),
"calibrated flux uncertainty",
"nJy"),
true);
312 newKey.magErr =
mapper.addOutputField(
313 FieldD(inSchema.join(
field,
"magErr"),
"calibrated magnitude uncertainty",
"mag(AB)"),
318 keys.emplace_back(newKey);
323 output.insert(
mapper, output.begin(), catalog.begin(), catalog.end());
325 auto calibration = evaluateCatalog(output);
329 for (
auto &rec : output) {
330 for (
auto &key :
keys) {
331 double instFlux = rec.get(key.instFlux);
332 double nanojansky = toNanojansky(instFlux, calibration[iRec]);
333 rec.set(key.flux, nanojansky);
334 rec.set(key.mag, toMagnitude(instFlux, calibration[iRec]));
335 if (key.instFluxErr.isValid()) {
336 double instFluxErr = rec.get(key.instFluxErr);
337 rec.set(key.fluxErr, toNanojanskyErr(instFlux, instFluxErr, calibration[iRec],
338 _calibrationErr, nanojansky));
340 toMagnitudeErr(instFlux, instFluxErr, calibration[iRec], _calibrationErr));
352 for (
auto const &
name : catalog.getSchema().getNames()) {
354 if (
name.size() > SUFFIX.
size() + 1 &&
359 return calibrateCatalog(catalog, instFluxFields);
366 class PhotoCalibSchema {
376 PhotoCalibSchema(PhotoCalibSchema
const &) =
delete;
377 PhotoCalibSchema &
operator=(PhotoCalibSchema
const &) =
delete;
379 PhotoCalibSchema(PhotoCalibSchema &&) =
delete;
380 PhotoCalibSchema &
operator=(PhotoCalibSchema &&) =
delete;
382 static PhotoCalibSchema
const &get() {
383 static PhotoCalibSchema
const instance;
391 "calibrationMean",
"mean calibration on this PhotoCalib's domain",
"count")),
393 schema.addField<double>(
"calibrationErr",
"1-sigma error on calibrationMean",
"count")),
394 isConstant(
schema.addField<table::Flag>(
"isConstant",
"Is this spatially-constant?")),
395 field(
schema.addField<int>(
"field",
"archive ID of the BoundedField object")),
396 version(
schema.addField<int>(
"version",
"version of this PhotoCalib")) {}
399 class PhotoCalibFactory :
public table::io::PersistableFactory {
402 read(InputArchive
const &archive, CatalogVector
const &catalogs)
const override {
403 table::BaseRecord
const &record = catalogs.front().front();
404 PhotoCalibSchema
const &
keys = PhotoCalibSchema::get();
405 int version = getVersion(record);
407 throw(pex::exceptions::RuntimeError(
"Unsupported version (version 0 was defined in maggies): " +
410 return std::make_shared<PhotoCalib>(record.get(
keys.calibrationMean), record.get(
keys.calibrationErr),
411 archive.get<afw::math::BoundedField>(record.get(
keys.field)),
412 record.get(
keys.isConstant));
418 int getVersion(table::BaseRecord
const &record)
const {
422 auto versionKey = record.getSchema().
find<
int>(versionName);
423 version = record.get(versionKey.key);
424 }
catch (
const pex::exceptions::NotFoundError &) {
432 std::string getPhotoCalibPersistenceName() {
return "PhotoCalib"; }
434 PhotoCalibFactory registration(getPhotoCalibPersistenceName());
443 int const CALIB_TABLE_CURRENT_VERSION = 2;
455 CalibKeys(
const CalibKeys &) =
delete;
456 CalibKeys &
operator=(
const CalibKeys &) =
delete;
459 CalibKeys(CalibKeys &&) =
delete;
460 CalibKeys &
operator=(CalibKeys &&) =
delete;
462 CalibKeys(
int tableVersion = CALIB_TABLE_CURRENT_VERSION)
464 if (tableVersion == 1) {
467 "midtime",
"middle of the time of the exposure relative to Unix epoch",
"ns");
475 class CalibFactory :
public table::io::PersistableFactory {
478 CatalogVector
const &catalogs)
const override {
481 int tableVersion = 1;
483 catalogs.front().getSchema().find<
double>(EXPTIME_FIELD_NAME);
484 }
catch (pex::exceptions::NotFoundError
const&) {
485 tableVersion = CALIB_TABLE_CURRENT_VERSION;
488 CalibKeys
const keys{tableVersion};
492 table::BaseRecord
const &record = catalogs.front().front();
503 std::string getCalibPersistenceName() {
return "Calib"; }
505 CalibFactory calibRegistration(getCalibPersistenceName());
512 PhotoCalibSchema
const &
keys = PhotoCalibSchema::get();
514 auto record = catalog.
addNew();
515 record->set(
keys.calibrationMean, _calibrationMean);
516 record->set(
keys.calibrationErr, _calibrationErr);
517 record->set(
keys.isConstant, _isConstant);
518 record->set(
keys.field, handle.
put(_calibration));
519 record->set(
keys.version, SERIALIZATION_VERSION);
527 return _calibrationMean;
529 return _calibration->evaluate(point);
532 ndarray::Array<double, 1> PhotoCalib::evaluateArray(ndarray::Array<double, 1>
const &xx,
533 ndarray::Array<double, 1>
const &yy)
const {
535 ndarray::Array<double, 1>
result = ndarray::allocate(ndarray::makeVector(xx.size()));
536 result.deep() = _calibrationMean;
539 return _calibration->evaluate(xx, yy);
544 ndarray::Array<double, 1> xx = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
545 ndarray::Array<double, 1> yy = ndarray::allocate(ndarray::makeVector(sourceCatalog.size()));
547 for (
auto const &rec : sourceCatalog) {
548 auto point = rec.getCentroid();
549 xx[i] = point.getX();
550 yy[i] = point.getY();
553 return evaluateArray(xx, yy);
558 ndarray::Array<double, 2, 2>
result)
const {
559 auto instFluxKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFlux").key;
560 auto instFluxErrKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFluxErr").key;
562 auto calibration = evaluateCatalog(sourceCatalog);
565 for (
auto const &rec : sourceCatalog) {
566 double instFlux = rec.get(instFluxKey);
567 double instFluxErr = rec.get(instFluxErrKey);
568 double nanojansky = toNanojansky(instFlux, calibration[i]);
569 (*iter)[0] = nanojansky;
570 (*iter)[1] = toNanojanskyErr(instFlux, instFluxErr, calibration[i], _calibrationErr, nanojansky);
578 ndarray::Array<double, 2, 2>
result)
const {
579 auto instFluxKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFlux").key;
580 auto instFluxErrKey = sourceCatalog.getSchema().find<
double>(instFluxField +
"_instFluxErr").key;
582 auto calibration = evaluateCatalog(sourceCatalog);
585 for (
auto const &rec : sourceCatalog) {
586 double instFlux = rec.get(instFluxKey);
587 double instFluxErr = rec.get(instFluxErrKey);
588 (*iter)[0] = toMagnitude(instFlux, calibration[i]);
589 (*iter)[1] = toMagnitudeErr(instFlux, instFluxErr, calibration[i], _calibrationErr);
596 auto key =
"FLUXMAG0";
597 if (metadata.
exists(key)) {
601 double instFluxMag0Err = 0.0;
603 if (metadata.
exists(key)) {
table::Key< std::string > name
#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
table::Key< int > version
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.
VariancePtr getVariance() const
Return a (shared_ptr to) the MaskedImage's variance.
ImagePtr getImage() const
Return a (shared_ptr to) the MaskedImage's image.
The photometric calibration of an exposure.
Schema getSchema() const
Return the Schema that holds this record's fields and keys.
Field< T >::Value get(Key< T > const &key) const
Return the value of a field for the given key.
std::shared_ptr< RecordT > addNew()
Create a new record, add it to the end of the catalog, and return a pointer to it.
Key< T > addField(Field< T > const &field, bool doReplace=false)
Add a new field to the Schema, and return the associated Key.
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).
Reports errors in the logical structure of the program.
Reports attempts to access elements using an invalid key.
T emplace_back(T... args)
def scale(algorithm, min, max=None, frame=None)
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.
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.
void write(OutputArchiveHandle &handle) const override
FilterProperty & operator=(FilterProperty const &)=default
std::ostream & operator<<(std::ostream &os, PhotoCalib const &photoCalib)
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.
Angle abs(Angle const &a)
A base class for image defects.
T setprecision(T... args)
A description of a field in a table.
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