Alexandria 2.25.0
SDC-CH common library for the Euclid project
PhotometryAttributeFromRow.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2012-2021 Euclid Science Ground Segment
3 *
4 * This library is free software; you can redistribute it and/or modify it under
5 * the terms of the GNU Lesser General Public License as published by the Free
6 * Software Foundation; either version 3.0 of the License, or (at your option)
7 * any later version.
8 *
9 * This library is distributed in the hope that it will be useful, but WITHOUT
10 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11 * FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
12 * details.
13 *
14 * You should have received a copy of the GNU Lesser General Public License
15 * along with this library; if not, write to the Free Software Foundation, Inc.,
16 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
17 */
18
27#include "ElementsKernel/Real.h"
28#include <cmath>
29#include <typeindex>
30
31#include "../../SourceCatalog/PhotometryParsingException.h"
33#include "ElementsKernel/Real.h"
34#include "Table/CastVisitor.h"
35
36namespace Euclid {
37namespace SourceCatalog {
38
42 const bool missing_photometry_enabled, const double missing_photometry_flag, const bool upper_limit_enabled,
43 const std::vector<std::pair<std::string, float>>& n_map, const double n_upper_limit_flag,
44 const std::vector<std::pair<std::string, bool>>& convert_from_mag)
45 : m_missing_photometry_enabled(missing_photometry_enabled)
46 , m_missing_photometry_flag(missing_photometry_flag)
47 , m_upper_limit_enabled(upper_limit_enabled)
48 , m_n_map(n_map)
49 , m_n_upper_limit_flag(n_upper_limit_flag)
50 , m_convert_from_mag(convert_from_mag) {
51
52 std::unique_ptr<size_t> flux_column_index_ptr;
53 std::unique_ptr<size_t> error_column_index_ptr;
54
55 for (auto filter_name_pair : filter_name_mapping) {
56 flux_column_index_ptr = column_info_ptr->find(filter_name_pair.second.first);
57 error_column_index_ptr = column_info_ptr->find(filter_name_pair.second.second);
58
59 if (flux_column_index_ptr == nullptr) {
60 throw Elements::Exception() << "Column info does not have the flux column " << filter_name_pair.second.first;
61 }
62 if (error_column_index_ptr == nullptr) {
63 throw Elements::Exception() << "Column info does not have the flux error column "
64 << filter_name_pair.second.second;
65 }
66 m_table_index_vector.push_back(std::make_pair(*(flux_column_index_ptr), *(error_column_index_ptr)));
67 }
68
69 // create and filled the shared pointer to the filter name vector
70 m_filter_name_vector_ptr = std::make_shared<std::vector<std::string>>();
71 for (auto a_filter_name_map : filter_name_mapping) {
72 m_filter_name_vector_ptr->push_back(a_filter_name_map.first);
73 }
74
75 // default the m_convert_from_mag
78 for (auto m_n_pair : m_n_map) {
79 m_convert_from_mag.push_back(std::make_pair(m_n_pair.first, false));
80 }
81 }
82}
83
85
86std::pair<double, double> PhotometryAttributeFromRow::convertFromMag(const double mag, const double mag_err) const {
87
90 } else {
91 // check if the error is a flag
92 bool is_flag = Elements::isEqual(mag_err, m_n_upper_limit_flag);
93
94 // compute the flux and the error
95 double flux = 3.631e9 * std::pow(10, -0.4 * mag);
96 // with this formula the sign of mag_err is forwarded to the flux_err
97 double flux_err = is_flag ? m_n_upper_limit_flag : 0.4 * flux * mag_err * std::log(10);
98
99 return std::make_pair(flux, flux_err);
100 }
101}
102
103static std::string getContextDescription(bool missing_photometry_enabled, bool upper_limit_enabled) {
104 std::string setup_desc("with ");
105 if (missing_photometry_enabled && upper_limit_enabled) {
106 setup_desc += "'missing data' and 'upper limit' enabled";
107 } else if (missing_photometry_enabled) {
108 setup_desc += "'missing data' and 'upper limit' disabled";
109 } else if (upper_limit_enabled) {
110 setup_desc += "'missing data' disabled and 'upper limit' enabled";
111 } else {
112 setup_desc += "'missing data' and 'upper limit' disabled";
113 }
114 return setup_desc;
115}
116
118 std::vector<FluxErrorPair> photometry_vector{};
119
121
122 auto n_threshod_iter = m_n_map.begin();
123 auto convert_from_mag_iter = m_convert_from_mag.begin();
124 for (auto& filter_index_pair : m_table_index_vector) {
125 Euclid::Table::Row::cell_type flux_cell = row[filter_index_pair.first];
126 Euclid::Table::Row::cell_type error_cell = row[filter_index_pair.second];
127
128 double flux = boost::apply_visitor(Table::CastVisitor<double>{}, flux_cell);
129 double error = boost::apply_visitor(Table::CastVisitor<double>{}, error_cell);
130
131 if (convert_from_mag_iter->second) {
132 auto converted = convertFromMag(flux, error);
133 flux = converted.first;
134 error = converted.second;
135 }
136
137 bool upper_limit = false;
138
139
140 bool missing_data = !std::isfinite(flux) | !std::isfinite(error) |
142
143 if (m_missing_photometry_enabled && missing_data) {
144 error = 0.;
145 } else if (!m_missing_photometry_enabled && missing_data) {
146 throw SourceCatalog::PhotometryParsingException("Non finite flux encountered when parsing the Photometry ",
147 context_desc.c_str(), flux, error);
148 } else if (m_upper_limit_enabled) {
150 if (error == 0.) {
151 throw SourceCatalog::PhotometryParsingException("Zero error encountered when parsing the Photometry",
152 context_desc.c_str(), flux, error);
153 }
154 if (error < 0) {
157 error = flux / n_threshod_iter->second;
158 }
159 upper_limit = true;
160 if (flux <= 0) {
162 "Negative or Zero flux encountered when parsing the Photometry", context_desc.c_str(), flux, error);
163 }
164
165 error = std::abs(error);
166 }
167 } else {
169 if (error <= 0) {
171 "Negative or Zero error encountered when parsing the Photometry", context_desc.c_str(), flux, error);
172 }
173 }
174
175 photometry_vector.push_back(FluxErrorPair{flux, error, missing_data, upper_limit});
176 ++n_threshod_iter;
177 ++convert_from_mag_iter;
178 } // Eof for
179
180 std::unique_ptr<Attribute> photometry_ptr{new Photometry{m_filter_name_vector_ptr, photometry_vector}};
181
182 return photometry_ptr;
183}
184
185} // namespace SourceCatalog
186} // end of namespace Euclid
T begin(T... args)
std::vector< std::pair< size_t, size_t > > m_table_index_vector
std::unique_ptr< Attribute > createAttribute(const Euclid::Table::Row &row) override
Create a photometricAttribute from a Table row.
std::pair< double, double > convertFromMag(const double mag, const double mag_err) const
std::vector< std::pair< std::string, bool > > m_convert_from_mag
std::shared_ptr< std::vector< std::string > > m_filter_name_vector_ptr
PhotometryAttributeFromRow(std::shared_ptr< Euclid::Table::ColumnInfo > column_info_ptr, const std::vector< std::pair< std::string, std::pair< std::string, std::string > > > &filter_name_mapping, const bool missing_photometry_enabled, const double missing_photometry_flag, const bool upper_limit_enabled, const std::vector< std::pair< std::string, float > > &n_map, const double n_upper_limit_flag, const std::vector< std::pair< std::string, bool > > &convert_from_mag={})
Create a PhotometryAttributeFromRow object.
std::vector< std::pair< std::string, float > > m_n_map
std::unique_ptr< std::size_t > find(const std::string &name) const
Returns the index of a column, given the name of it, or nullptr if there is no column with this name.
Definition: ColumnInfo.cpp:77
Represents one row of a Table.
Definition: Row.h:64
boost::variant< bool, int32_t, int64_t, float, double, std::string, std::vector< bool >, std::vector< int32_t >, std::vector< int64_t >, std::vector< float >, std::vector< double >, NdArray::NdArray< int32_t >, NdArray::NdArray< int64_t >, NdArray::NdArray< float >, NdArray::NdArray< double > > cell_type
The possible cell types.
Definition: Row.h:71
T isfinite(T... args)
T isnan(T... args)
T log(T... args)
T make_pair(T... args)
bool isEqual(const RawType &left, const RawType &right)
static std::string getContextDescription(bool missing_photometry_enabled, bool upper_limit_enabled)
T pow(T... args)
T push_back(T... args)
T size(T... args)