Rolf Buehler / gammalib

Commit e8e73b7743f2a06ef2558d78dc3dfcb532774972

Authored by Jürgen Knödlseder
1 parent 8b6160d2

Handle appending of additional columns to GArf structure.

Inline Side-by-side
Showing 6 changed files with 295 additions and 13 deletions
ChangeLog
  View file @e8e73b7
... ... @@ -3,6 +3,7 @@
3 3 * Version 1.1.0 released
4 4 =========================
5 5  
  6 + Handle additional columns in GArf class
6 7 GCTAAeff2D no longer extrapolates beyond response table (#1784)
7 8 Add GWcsSIN projection class (#1798)
8 9 CTA cube response no longer selects on RoI and energies (#1789)
... ...
... ... @@ -38,6 +38,11 @@ The following methods have been added to the public interface:
38 38 - GModelSky::spatial()
39 39 - GModelSky::spectral()
40 40 - GModelSky::temporal()
  41 +- GArf::columns()
  42 +- GArf::at()
  43 +- GArf::append()
  44 +- GArf::operator()(int,int)
  45 +- GArf::operator[](std::string)
41 46 - GCTAObservation(GFilename) constructor
42 47 - GCTAResponseTable::axis(std::string)
43 48 - GCTAResponseTable::table(std::string)
... ... @@ -308,7 +313,8 @@ location and the Hub will be informed that an image is available.
308 313  
309 314 16. Xspec module
310 315 ----------------
311   -None
  316 +The GArf class now also handles additional double precision columns that may
  317 +be needed for response computation.
312 318  
313 319  
314 320 17. Test module
... ...
include/GArf.hpp
  View file @e8e73b7
... ... @@ -58,14 +58,25 @@ public:
58 58 GArf& operator=(const GArf& arf);
59 59 double& operator[](const int& index);
60 60 const double& operator[](const int& index) const;
  61 + double& operator()(const int& index, const int& col);
  62 + const double& operator()(const int& index, const int& col) const;
  63 +
  64 + // Additional column access operators
  65 + std::vector<double>& operator[](const std::string& colname);
  66 + const std::vector<double>& operator[](const std::string& colname) const;
61 67  
62 68 // Methods
63 69 void clear(void);
64 70 GArf* clone(void) const;
65 71 std::string classname(void) const;
66 72 int size(void) const;
  73 + int columns(void) const;
67 74 double& at(const int& index);
68 75 const double& at(const int& index) const;
  76 + double& at(const int& index, const int& col);
  77 + const double& at(const int& index, const int& col) const;
  78 + void append(const std::string& name,
  79 + const std::vector<double>& column);
69 80 const GEbounds& ebounds(void) const;
70 81 void load(const GFilename& filename);
71 82 void save(const GFilename& filename,
... ... @@ -80,11 +91,14 @@ protected:
80 91 void init_members(void);
81 92 void copy_members(const GArf& pha);
82 93 void free_members(void);
  94 + int column_index(const std::string& colname) const;
83 95  
84 96 // Protected members
85   - mutable GFilename m_filename; //!< Filename of origin
86   - GEbounds m_ebounds; //!< Energy boundaries
87   - std::vector<double> m_specresp; //!< Spectral response
  97 + mutable GFilename m_filename; //!< Filename of origin
  98 + GEbounds m_ebounds; //!< Energy boundaries
  99 + std::vector<double> m_specresp; //!< Spectral response
  100 + std::vector<std::string> m_colnames; //!< Additional column names
  101 + std::vector<std::vector<double> > m_coldata; //!< Additional column data
88 102 };
89 103  
90 104  
... ... @@ -129,6 +143,36 @@ const double&amp; GArf::operator[](const int&amp; index) const
129 143  
130 144  
131 145 /***********************************************************************//**
  146 + * @brief Return content of additional columns
  147 + *
  148 + * @param[in] index Bin index [0,...,size()-1].
  149 + * @param[in] col Columns index [0,...,columns()-1].
  150 + *
  151 + * Returns reference to content of additional columns.
  152 + ***************************************************************************/
  153 +inline
  154 +double& GArf::operator()(const int& index, const int& col)
  155 +{
  156 + return (m_coldata[col][index]);
  157 +}
  158 +
  159 +
  160 +/***********************************************************************//**
  161 + * @brief Return content of additional columns (const version)
  162 + *
  163 + * @param[in] index Bin index [0,...,size()-1].
  164 + * @param[in] col Columns index [0,...,columns()-1].
  165 + *
  166 + * Returns reference to content of additional columns.
  167 + ***************************************************************************/
  168 +inline
  169 +const double& GArf::operator()(const int& index, const int& col) const
  170 +{
  171 + return (m_coldata[col][index]);
  172 +}
  173 +
  174 +
  175 +/***********************************************************************//**
132 176 * @brief Return number of spectral bins
133 177 *
134 178 * @return Number of spectral bins.
... ... @@ -143,6 +187,20 @@ int GArf::size(void) const
143 187  
144 188  
145 189 /***********************************************************************//**
  190 + * @brief Return number of additional columns
  191 + *
  192 + * @return Number of additional columns.
  193 + *
  194 + * Returns the number of additional columns.
  195 + ***************************************************************************/
  196 +inline
  197 +int GArf::columns(void) const
  198 +{
  199 + return (int)m_colnames.size();
  200 +}
  201 +
  202 +
  203 +/***********************************************************************//**
146 204 * @brief Return energy boundaries
147 205 *
148 206 * @return Energy boundaries for all spectral bins.
... ...
pyext/GArf.i
  View file @e8e73b7
... ... @@ -44,12 +44,20 @@ public:
44 44 GArf(const GArf& arf);
45 45 virtual ~GArf(void);
46 46  
  47 + // Operators
  48 + double& operator()(const int& index, const int& col);
  49 + const double& operator()(const int& index, const int& col) const;
  50 +
47 51 // Methods
48 52 void clear(void);
49 53 GArf* clone(void) const;
50 54 std::string classname(void) const;
51 55 int size(void) const;
  56 + int columns(void) const;
52 57 double& at(const int& index);
  58 + double& at(const int& index, const int& col);
  59 + void append(const std::string& name,
  60 + const std::vector<double>& column);
53 61 const GEbounds& ebounds(void) const;
54 62 void load(const GFilename& filename);
55 63 void save(const GFilename& filename,
... ...
src/xspec/GArf.cpp
  View file @e8e73b7
... ... @@ -39,7 +39,10 @@
39 39 #include "GFitsTableFloatCol.hpp"
40 40  
41 41 /* __ Method name definitions ____________________________________________ */
42   -#define G_AT "GArf::at(int&)"
  42 +#define G_OPERATOR "GArf::operator[](std::string&)"
  43 +#define G_AT1 "GArf::at(int&)"
  44 +#define G_AT2 "GArf::at(int&, int&)"
  45 +#define G_APPEND "GArf::append(std::string&, std::vector<double>&)"
43 46  
44 47 /* __ Macros _____________________________________________________________ */
45 48  
... ... @@ -170,6 +173,56 @@ GArf&amp; GArf::operator=(const GArf&amp; arf)
170 173 }
171 174  
172 175  
  176 +/***********************************************************************//**
  177 + * @brief Return additional vector column
  178 + *
  179 + * @param[in] colname Vector column name.
  180 + * @return Vector column.
  181 + *
  182 + * Returns reference to additional vector column.
  183 + ***************************************************************************/
  184 +std::vector<double>& GArf::operator[](const std::string& colname)
  185 +{
  186 + // Determine index of additional column (-1 if not found)
  187 + int index = column_index(colname);
  188 +
  189 + // Throw exception if index not found
  190 + if (index == -1) {
  191 + std::string msg = "Could not find additional column with name \""+
  192 + colname+"\".";
  193 + throw GException::invalid_value(G_OPERATOR, msg);
  194 + }
  195 +
  196 + // Return reference to vector column
  197 + return (m_coldata[index]);
  198 +}
  199 +
  200 +
  201 +/***********************************************************************//**
  202 + * @brief Return additional vector column (const version)
  203 + *
  204 + * @param[in] colname Vector column name.
  205 + * @return Vector column.
  206 + *
  207 + * Returns reference to additional vector column.
  208 + ***************************************************************************/
  209 +const std::vector<double>& GArf::operator[](const std::string& colname) const
  210 +{
  211 + // Determine index of additional column (-1 if not found)
  212 + int index = column_index(colname);
  213 +
  214 + // Throw exception if index not found
  215 + if (index == -1) {
  216 + std::string msg = "Could not find additional column with name \""+
  217 + colname+"\".";
  218 + throw GException::invalid_value(G_OPERATOR, msg);
  219 + }
  220 +
  221 + // Return reference to vector column
  222 + return (m_coldata[index]);
  223 +}
  224 +
  225 +
173 226 /*==========================================================================
174 227 = =
175 228 = Public methods =
... ... @@ -216,9 +269,9 @@ GArf* GArf::clone(void) const
216 269 ***************************************************************************/
217 270 double& GArf::at(const int& index)
218 271 {
219   - // Raise exception if index is out of range
  272 + // Throw an exception if index is out of range
220 273 if (index < 0 || index >= size()) {
221   - throw GException::out_of_range(G_AT, index, 0, size()-1);
  274 + throw GException::out_of_range(G_AT1, index, 0, size()-1);
222 275 }
223 276  
224 277 // Return reference
... ... @@ -238,9 +291,9 @@ double&amp; GArf::at(const int&amp; index)
238 291 ***************************************************************************/
239 292 const double& GArf::at(const int& index) const
240 293 {
241   - // Raise exception if index is out of range
  294 + // Throw an exception if index is out of range
242 295 if (index < 0 || index >= size()) {
243   - throw GException::out_of_range(G_AT, index, 0, size()-1);
  296 + throw GException::out_of_range(G_AT1, index, 0, size()-1);
244 297 }
245 298  
246 299 // Return reference
... ... @@ -249,6 +302,86 @@ const double&amp; GArf::at(const int&amp; index) const
249 302  
250 303  
251 304 /***********************************************************************//**
  305 + * @brief Return content of additional columns
  306 + *
  307 + * @param[in] index Bin index [0,...,size()-1].
  308 + * @param[in] col Columns index [0,...,columns()-1].
  309 + *
  310 + * @exception GException::out_of_range
  311 + * Spectral bin or column index is out of range.
  312 + *
  313 + * Returns reference to content of additional columns.
  314 + ***************************************************************************/
  315 +double& GArf::at(const int& index, const int& col)
  316 +{
  317 + // Throw an exception if bin or column index is out of range
  318 + if (index < 0 || index >= size()) {
  319 + throw GException::out_of_range(G_AT2, "Bin index", index, size());
  320 + }
  321 + if (col < 0 || col >= columns()) {
  322 + throw GException::out_of_range(G_AT2, "Column index", col, columns());
  323 + }
  324 +
  325 + // Return reference
  326 + return (m_coldata[col][index]);
  327 +}
  328 +
  329 +
  330 +/***********************************************************************//**
  331 + * @brief Return content of additional columns (const version)
  332 + *
  333 + * @param[in] index Bin index [0,...,size()-1].
  334 + * @param[in] col Columns index [0,...,columns()-1].
  335 + *
  336 + * @exception GException::out_of_range
  337 + * Spectral bin or column index is out of range.
  338 + *
  339 + * Returns reference to content of additional columns.
  340 + ***************************************************************************/
  341 +const double& GArf::at(const int& index, const int& col) const
  342 +{
  343 + // Throw an exception if bin or column index is out of range
  344 + if (index < 0 || index >= size()) {
  345 + throw GException::out_of_range(G_AT2, "Bin index", index, size());
  346 + }
  347 + if (col < 0 || col >= columns()) {
  348 + throw GException::out_of_range(G_AT2, "Column index", col, columns());
  349 + }
  350 +
  351 + // Return reference
  352 + return (m_coldata[col][index]);
  353 +}
  354 +
  355 +
  356 +/***********************************************************************//**
  357 + * @brief Append additional column to spectrum
  358 + *
  359 + * @param[in] name Additional column name.
  360 + * @param[in] column Additional column data.
  361 + ***************************************************************************/
  362 +void GArf::append(const std::string& name, const std::vector<double>& column)
  363 +{
  364 + // Throw an exception if the number of elements in the column does not
  365 + // correspond to the size of the spectrum
  366 + if (column.size() != size()) {
  367 + std::string msg = "Size of column "+gammalib::str(column.size())+
  368 + " is incompatible with size of spectrum "+
  369 + gammalib::str(size())+".";
  370 + throw GException::invalid_argument(G_APPEND, msg);
  371 + }
  372 +
  373 + // Append column name
  374 + m_colnames.push_back(name);
  375 +
  376 + // Append column data
  377 + m_coldata.push_back(column);
  378 +
  379 + // Return
  380 + return;
  381 +}
  382 +
  383 +
  384 +/***********************************************************************//**
252 385 * @brief Load Auxiliary Response File
253 386 *
254 387 * @param[in] filename File name.
... ... @@ -346,7 +479,8 @@ void GArf::read(const GFitsTable&amp; table)
346 479  
347 480 // Determine effective area conversion factor. Internal
348 481 // units are cm^2
349   - std::string u_specresp = gammalib::tolower(gammalib::strip_whitespace(specresp->unit()));
  482 + std::string u_specresp =
  483 + gammalib::tolower(gammalib::strip_whitespace(specresp->unit()));
350 484 double c_specresp = 1.0;
351 485 if (u_specresp == "m^2" || u_specresp == "m2") {
352 486 c_specresp = 10000.0;
... ... @@ -369,6 +503,31 @@ void GArf::read(const GFitsTable&amp; table)
369 503  
370 504 } // endfor: looped over energy bins
371 505  
  506 + // Read any additional columns
  507 + for (int icol = 0; icol < table.ncols(); ++icol) {
  508 +
  509 + // Fall through if the column is a standard column
  510 + std::string colname(table[icol]->name());
  511 + if ((colname == "ENERG_LO") ||
  512 + (colname == "ENERG_HI") ||
  513 + (colname == "SPECRESP")) {
  514 + continue;
  515 + }
  516 +
  517 + // Get pointer to column
  518 + const GFitsTableCol* column = table[icol];
  519 +
  520 + // Set column vector
  521 + std::vector<double> coldata;
  522 + for (int i = 0; i < num; ++i) {
  523 + coldata.push_back(column->real(i));
  524 + }
  525 +
  526 + // Append column
  527 + append(colname, coldata);
  528 +
  529 + } // endfor: looped over all additional columns
  530 +
372 531 // Return
373 532 return;
374 533 }
... ... @@ -415,9 +574,9 @@ void GArf::write(GFits&amp; fits) const
415 574  
416 575 // Fill columns
417 576 for (int i = 0; i < length; ++i) {
418   - energy_lo(i) = m_ebounds.emin(i).keV();
419   - energy_hi(i) = m_ebounds.emax(i).keV();
420   - specresp(i) = m_specresp[i];
  577 + energy_lo(i) = (float)m_ebounds.emin(i).keV();
  578 + energy_hi(i) = (float)m_ebounds.emax(i).keV();
  579 + specresp(i) = (float)m_specresp[i];
421 580 }
422 581  
423 582 // Set column units
... ... @@ -433,6 +592,22 @@ void GArf::write(GFits&amp; fits) const
433 592 hdu.append(energy_hi);
434 593 hdu.append(specresp);
435 594  
  595 + // Append any additional columns
  596 + for (int icol = 0; icol < columns(); ++icol) {
  597 +
  598 + // Allocate floating point vector columns
  599 + GFitsTableFloatCol column(m_colnames[icol], length);
  600 +
  601 + // Fill columns
  602 + for (int i = 0; i < length; ++i) {
  603 + column(i) = (float)m_coldata[icol][i];
  604 + }
  605 +
  606 + // Append columns to table
  607 + hdu.append(column);
  608 +
  609 + } // endfor: looped over all additional columns
  610 +
436 611 // Append HDU to FITS file
437 612 fits.append(hdu);
438 613  
... ... @@ -490,6 +665,8 @@ void GArf::init_members(void)
490 665 m_filename.clear();
491 666 m_ebounds.clear();
492 667 m_specresp.clear();
  668 + m_colnames.clear();
  669 + m_coldata.clear();
493 670  
494 671 // Return
495 672 return;
... ... @@ -507,6 +684,8 @@ void GArf::copy_members(const GArf&amp; arf)
507 684 m_filename = arf.m_filename;
508 685 m_ebounds = arf.m_ebounds;
509 686 m_specresp = arf.m_specresp;
  687 + m_colnames = arf.m_colnames;
  688 + m_coldata = arf.m_coldata;
510 689  
511 690 // Return
512 691 return;
... ... @@ -521,3 +700,30 @@ void GArf::free_members(void)
521 700 // Return
522 701 return;
523 702 }
  703 +
  704 +
  705 +/***********************************************************************//**
  706 + * @brief Returns index of additional vector column
  707 + *
  708 + * @param[in] colname Vector column name.
  709 + * @return Vector column index.
  710 + *
  711 + * Returns index of additional vector column. Returns -1 if the vector column
  712 + * has not found.
  713 + ***************************************************************************/
  714 +int GArf::column_index(const std::string& colname) const
  715 +{
  716 + // Initialise index
  717 + int index(-1);
  718 +
  719 + // Search vector column name
  720 + for (int i = 0; i < columns(); ++i) {
  721 + if (m_colnames[i] == colname) {
  722 + index = i;
  723 + break;
  724 + }
  725 + }
  726 +
  727 + // Return index
  728 + return index;
  729 +}
... ...
src/xspec/GPha.cpp
  View file @e8e73b7
... ... @@ -425,6 +425,9 @@ void GPha::save(const GFilename&amp; filename, const bool&amp; clobber) const
425 425 ***************************************************************************/
426 426 void GPha::read(const GFitsTable& table)
427 427 {
  428 + // Clear spectrum
  429 + clear();
  430 +
428 431 // Get data column
429 432 const GFitsTableCol* col_data = table["COUNTS"];
430 433  
... ...