GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAResponseTable.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAResponseTable.cpp - CTA response table class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2018 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCTAResponseTable.cpp
23  * @brief CTA response table class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GException.hpp"
35 #include "GFitsBinTable.hpp"
36 #include "GFitsTableFloatCol.hpp"
37 #include "GCTAResponseTable.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_OPERATOR1 "GCTAResponseTable::operator()(double&)"
41 #define G_OPERATOR2 "GCTAResponseTable::operator()(double&, double&)"
42 #define G_OPERATOR3 "GCTAResponseTable::operator()(double&, double&,"\
43  " double&)"
44 #define G_ELEMENT_OPERATOR1 "GCTAResponseTable::operator()(int&)"
45 #define G_ELEMENT_OPERATOR2 "GCTAResponseTable::operator()(int&, int&)"
46 #define G_INX_OPERATOR1 "GCTAResponseTable::operator()(int&, double&)"
47 #define G_INX_OPERATOR2 "GCTAResponseTable::operator()(int&, double&,"\
48  " double&)"
49 #define G_INX_OPERATOR3 "GCTAResponseTable::operator()(int&, double&,"\
50  " double&, double&)"
51 #define G_TABLE "GCTAResponseTable::table(std::string&)"
52 #define G_SCALE "GCTAResponseTable::scale(int&, double&)"
53 #define G_AXIS "GCTAResponseTable::axis(std::string&)"
54 #define G_AXIS_BINS "GCTAResponseTable::axis_bins(int&)"
55 #define G_AXIS_LO_NAME "GCTAResponseTable::axis_lo_name(int&)"
56 #define G_AXIS_HI_NAME "GCTAResponseTable::axis_hi_name(int&)"
57 #define G_AXIS_LO_UNIT "GCTAResponseTable::axis_lo_unit(int&)"
58 #define G_AXIS_HI_UNIT "GCTAResponseTable::axis_hi_unit(int&)"
59 #define G_UNIT "GCTAResponseTable::unit(int&)"
60 #define G_AXIS_LO "GCTAResponseTable::axis_lo(int&, int&)"
61 #define G_AXIS_HI "GCTAResponseTable::axis_hi(int&, int&)"
62 #define G_AXIS_LINEAR "GCTAResponseTable::axis_linear(int&)"
63 #define G_AXIS_LOG10 "GCTAResponseTable::axis_log10(int&)"
64 #define G_AXIS_RADIANS "GCTAResponseTable::axis_radians(int&)"
65 #define G_AXIS_NODES "GCTAResponseTable::axis_nodes(int&)"
66 #define G_APPEND_AXIS "GCTAResponseTable::append_axis(std::vector<double>&,"\
67  " std::vector<double>&, std::string&, std::string&)"
68 #define G_APPEND_TABLE "GCTAResponseTable::append_table(std::string&,"\
69  " std::string&)"
70 #define G_READ "GCTAResponseTable::read(GFitsTable&)"
71 #define G_READ_COLNAMES "GCTAResponseTable::read_colnames(GFitsTable&)"
72 #define G_READ_AXES "GCTAResponseTable::read_axes(GFitsTable&)"
73 #define G_READ_TABLES "GCTAResponseTable::read_tables(GFitsTable&)"
74 
75 /* __ Macros _____________________________________________________________ */
76 
77 /* __ Coding definitions _________________________________________________ */
78 
79 /* __ Debug definitions __________________________________________________ */
80 
81 /* __ Constants __________________________________________________________ */
82 
83 
84 /*==========================================================================
85  = =
86  = Constructors/destructors =
87  = =
88  ==========================================================================*/
89 
90 /***********************************************************************//**
91  * @brief Void constructor
92  *
93  * Construct an empty CTA response table.
94  ***************************************************************************/
96 {
97  // Initialise class members
98  init_members();
99 
100  // Return
101  return;
102 }
103 
104 
105 /***********************************************************************//**
106  * @brief Copy constructor
107  *
108  * @param[in] table Response table.
109  *
110  * Construct a CTA response table by copying information from an existing
111  * response table. A deep copy is performed for construction, so that the
112  * original object may be destroyed after construction the new object.
113  ***************************************************************************/
115 {
116  // Initialise class members
117  init_members();
118 
119  // Copy members
120  copy_members(table);
121 
122  // Return
123  return;
124 }
125 
126 
127 /***********************************************************************//**
128  * @brief FITS table constructor
129  *
130  * @param[in] hdu FITS table HDU.
131  *
132  * Construct a CTA response table from the information in a FITS table HDU.
133  * Refer to the read() method for details.
134  ***************************************************************************/
136 {
137  // Initialise class members
138  init_members();
139 
140  // Read table information from HDU
141  read(hdu);
142 
143  // Return
144  return;
145 }
146 
147 
148 /***********************************************************************//**
149  * @brief Destructor
150  *
151  * Destroys the CTA response table.
152  ***************************************************************************/
154 {
155  // Free members
156  free_members();
157 
158  // Return
159  return;
160 }
161 
162 
163 /*==========================================================================
164  = =
165  = Operators =
166  = =
167  ==========================================================================*/
168 
169 /***********************************************************************//**
170  * @brief Assignment operator
171  *
172  * @param[in] table Response table.
173  * @return Response table.
174  *
175  * Assigns a CTA response table to another object. The method is assigning
176  * a deep copy of the response table, so that the original object can be
177  * destroyed after assignment without any loss of information.
178  ***************************************************************************/
180 {
181  // Execute only if object is not identical
182  if (this != &table) {
183 
184  // Free members
185  free_members();
186 
187  // Initialise private members
188  init_members();
189 
190  // Copy members
191  copy_members(table);
192 
193  } // endif: object was not identical
194 
195  // Return this object
196  return *this;
197 }
198 
199 
200 /***********************************************************************//**
201  * @brief Linear interpolation operator for 1D tables
202  *
203  * @param[in] arg Value.
204  * @return Vector of linearly interpolated response tables.
205  *
206  * @exception GException::invalid_value
207  * Invalid response table dimension.
208  *
209  * Evaluates all response tables at a given value. This method applies to
210  * one-dimensional tables only. If the tables have a different dimension
211  * an exception is thrown.
212  *
213  * The evaluation is performed by a linear interpolation of the table. If
214  * the specified value lies outside the range covered by the table, the
215  * table is linearly extrapolated from using either the first or the last
216  * two vector elements.
217  *
218  * @todo Write down formula.
219  ***************************************************************************/
220 std::vector<double> GCTAResponseTable::operator()(const double& arg) const
221 {
222  // Throw exception if table is not 1D
223  if (axes() != 1) {
224  std::string msg = "Invalid response table dimension "+
225  gammalib::str(axes())+" encountered. Response "
226  "table needs to be one-dimensional.";
228  }
229 
230  // Initialise result vector
231  std::vector<double> result(m_ntables);
232 
233  // Set indices and weighting factors for interpolation
234  update(arg);
235 
236  // Perform 1D interpolation
237  for (int i = 0; i < m_ntables; ++i) {
238  result[i] = m_wgt_left * m_tables[i][m_inx_left] +
240  }
241 
242  // Return result vector
243  return result;
244 }
245 
246 
247 /***********************************************************************//**
248  * @brief Bilinear interpolation operator for 2D tables
249  *
250  * @param[in] arg1 Value for first axis.
251  * @param[in] arg2 Value for second axis.
252  * @return Vector of bilinearly interpolated response tables.
253  *
254  * @exception GException::invalid_value
255  * Invalid response table dimension.
256  *
257  * Evaluates all response tables at a given value. This method applies to
258  * two-dimensional tables only. If the tables have a different dimension
259  * an exception is thrown.
260  *
261  * The evaluation is performed by a bilinear interpolation of the table. If
262  * the specified value lies outside the range covered by the table, the
263  * table is linearly extrapolated from using either the first or the last
264  * two vector elements.
265  *
266  * @todo Write down formula.
267  ***************************************************************************/
268 std::vector<double> GCTAResponseTable::operator()(const double& arg1,
269  const double& arg2) const
270 {
271  // Throw exception if table is not 2D
272  if (axes() != 2) {
273  std::string msg = "Invalid response table dimension "+
274  gammalib::str(axes())+" encountered. Response "
275  "table needs to be two-dimensional.";
277  }
278 
279  // Initialise result vector
280  std::vector<double> result(m_ntables);
281 
282  // Set indices and weighting factors for interpolation
283  update(arg1, arg2);
284 
285  // Perform 2D interpolation
286  for (int i = 0; i < m_ntables; ++i) {
287  result[i] = m_wgt1 * m_tables[i][m_inx1] +
288  m_wgt2 * m_tables[i][m_inx2] +
289  m_wgt3 * m_tables[i][m_inx3] +
290  m_wgt4 * m_tables[i][m_inx4];
291  }
292 
293  // Return result vector
294  return result;
295 }
296 
297 
298 /***********************************************************************//**
299  * @brief Trilinear interpolation operator for 3D tables
300  *
301  * @param[in] arg1 Value for first axis.
302  * @param[in] arg2 Value for second axis.
303  * @param[in] arg3 Value for third axis.
304  * @return Vector of trilinearly interpolated response tables.
305  *
306  * @exception GException::invalid_value
307  * Invalid response table dimension.
308  *
309  * Evaluates all response tables at a given value. This method applies to
310  * three-dimensional tables only. If the tables have a different dimension
311  * an exception is thrown.
312  *
313  * The evaluation is performed by a trilinear interpolation of the table. If
314  * the specified value lies outside the range covered by the table, the
315  * table is linearly extrapolated from using either the first or the last
316  * two vector elements.
317  *
318  * @todo Write down formula.
319  ***************************************************************************/
320 std::vector<double> GCTAResponseTable::operator()(const double& arg1,
321  const double& arg2,
322  const double& arg3) const
323 {
324  // Throw exception if table is not 3D
325  if (axes() != 3) {
326  std::string msg = "Invalid response table dimension "+
327  gammalib::str(axes())+" encountered. Response "
328  "table needs to be tri-dimensional.";
330  }
331 
332  // Initialise result vector
333  std::vector<double> result(m_ntables);
334 
335  // Set indices and weighting factors for interpolation
336  update(arg1, arg2, arg3);
337 
338  // Perform 3D interpolation
339  for (int i = 0; i < m_ntables; ++i) {
340  result[i] = m_wgt1 * m_tables[i][m_inx1] +
341  m_wgt2 * m_tables[i][m_inx2] +
342  m_wgt3 * m_tables[i][m_inx3] +
343  m_wgt4 * m_tables[i][m_inx4] +
344  m_wgt5 * m_tables[i][m_inx5] +
345  m_wgt6 * m_tables[i][m_inx6] +
346  m_wgt7 * m_tables[i][m_inx7] +
347  m_wgt8 * m_tables[i][m_inx8] ;
348  }
349 
350  // Return result vector
351  return result;
352 }
353 
354 
355 /***********************************************************************//**
356  * @brief Table element access operator
357  *
358  * @param[in] element Table element index [0,...,elements()-1].
359  * @return Table element value.
360  *
361  * @exception GException::out_of_range
362  * Invalid element index.
363  *
364  * Returns the value of a table element for the first response table.
365  ***************************************************************************/
366 double& GCTAResponseTable::operator()(const int& element)
367 {
368  // Optionally check if the index and element is valid
369  #if defined(G_RANGE_CHECK)
370  if (element < 0 || element >= elements()) {
371  throw GException::out_of_range(G_ELEMENT_OPERATOR1, "Element index",
372  element, elements());
373  }
374  #endif
375 
376  // Return elements
377  return (m_tables[0][element]);
378 }
379 
380 
381 /***********************************************************************//**
382  * @brief Table element access operator (const version)
383  *
384  * @param[in] element Table element index [0,...,elements()-1].
385  * @return Table element value.
386  *
387  * @exception GException::out_of_range
388  * Invalid element index.
389  *
390  * Returns the value of a table element for the first response table.
391  ***************************************************************************/
392 const double& GCTAResponseTable::operator()(const int& element) const
393 {
394  // Optionally check if the index and element is valid
395  #if defined(G_RANGE_CHECK)
396  if (element < 0 || element >= elements()) {
397  throw GException::out_of_range(G_ELEMENT_OPERATOR1, "Element index",
398  element, elements());
399  }
400  #endif
401 
402  // Return elements
403  return (m_tables[0][element]);
404 }
405 
406 
407 /***********************************************************************//**
408  * @brief Table element access operator
409  *
410  * @param[in] table Table index [0,...,tables()-1].
411  * @param[in] element Element index [0,...,elements()-1].
412  * @return Response table element.
413  *
414  * @exception GException::out_of_range
415  * Invalid table or element index.
416  *
417  * Returns the value of a table element for the response table with index
418  * @p table.
419  ***************************************************************************/
420 double& GCTAResponseTable::operator()(const int& table, const int& element)
421 {
422  // Optionally check if the index and element is valid
423  #if defined(G_RANGE_CHECK)
424  if (table < 0 || table >= tables()) {
425  throw GException::out_of_range(G_ELEMENT_OPERATOR2, "Table index",
426  table, tables());
427  }
428  if (element < 0 || element >= elements()) {
429  throw GException::out_of_range(G_ELEMENT_OPERATOR2, "Element index",
430  element, elements());
431  }
432  #endif
433 
434  // Return elements
435  return (m_tables[table][element]);
436 }
437 
438 
439 /***********************************************************************//**
440  * @brief Table element access operator (const version)
441  *
442  * @param[in] table Table index [0,...,tables()-1].
443  * @param[in] element Element index [0,...,elements()-1].
444  * @return Response table element.
445  *
446  * @exception GException::out_of_range
447  * Invalid table or element index.
448  *
449  * Returns the value of a table element for the response table with index
450  * @p table.
451  ***************************************************************************/
452 const double& GCTAResponseTable::operator()(const int& table, const int& element) const
453 {
454  // Optionally check if the index and element is valid
455  #if defined(G_RANGE_CHECK)
456  if (table < 0 || table >= tables()) {
457  throw GException::out_of_range(G_ELEMENT_OPERATOR2, "Table index",
458  table, tables());
459  }
460  if (element < 0 || element >= elements()) {
461  throw GException::out_of_range(G_ELEMENT_OPERATOR2, "Element index",
462  element, elements());
463  }
464  #endif
465 
466  // Return elements
467  return (m_tables[table][element]);
468 }
469 
470 
471 /***********************************************************************//**
472  * @brief Linear interpolation operator for 1D tables
473  *
474  * @param[in] table Table index [0,...,tables()-1].
475  * @param[in] arg Value.
476  * @return Linearly interpolated response table.
477  *
478  * @exception GException::invalid_value
479  * Invalid response table dimension.
480  * @exception GException::out_of_range
481  * Invalid table index.
482  *
483  * Evaluates one response table at a given value. This method applies to
484  * one-dimensional response tables. If the table has a different dimension
485  * an exception is thrown.
486  *
487  * The evaluation is performed by a linear interpolation of the table values.
488  * If the specified value lies outside the range covered by the table, the
489  * table is linearly extrapolated from using either the first or the last
490  * two vector elements.
491  *
492  * @todo Write down formula.
493  ***************************************************************************/
494 double GCTAResponseTable::operator()(const int& table, const double& arg) const
495 {
496  // Throw exception if table is not 1D
497  if (axes() != 1) {
498  std::string msg = "Invalid response table dimension "+
499  gammalib::str(axes())+" encountered. Response "
500  "table needs to be one-dimensional.";
502  }
503 
504  // Throw exception if index is not valid
505  #if defined(G_RANGE_CHECK)
506  if (table < 0 || table >= tables()) {
507  throw GException::out_of_range(G_INX_OPERATOR1, "Table index",
508  table, tables());
509  }
510  #endif
511 
512  // Set indices and weighting factors for interpolation
513  update(arg);
514 
515  // Perform 1D interpolation
516  double result = m_wgt_left * m_tables[table][m_inx_left] +
518 
519  // Return result
520  return result;
521 }
522 
523 
524 /***********************************************************************//**
525  * @brief Bilinear interpolation operator for 2D tables
526  *
527  * @param[in] table Table index [0,...,tables()-1].
528  * @param[in] arg1 Value for first axis.
529  * @param[in] arg2 Value for second axis.
530  * @return Bilinearly interpolated response table.
531  *
532  * @exception GException::invalid_value
533  * Invalid response table dimension.
534  * @exception GException::out_of_range
535  * Invalid table index.
536  *
537  * Evaluates one response table at a given value. This method applies to
538  * two-dimensional response tables. If the table has a different dimension
539  * an exception is thrown.
540  *
541  * The evaluation is performed by a bilinear interpolation of the table
542  * values. If the specified value lies outside the range covered by the
543  * table, the table is linearly extrapolated from using either the first
544  * or the last two vector elements.
545  *
546  * @todo Write down formula.
547  ***************************************************************************/
548 double GCTAResponseTable::operator()(const int& table,
549  const double& arg1,
550  const double& arg2) const
551 {
552  // Throw exception if table is not 2D
553  if (axes() != 2) {
554  std::string msg = "Invalid response table dimension "+
555  gammalib::str(axes())+" encountered. Response "
556  "table needs to be two-dimensional.";
558  }
559 
560  // Throw exception if index is not valid
561  #if defined(G_RANGE_CHECK)
562  if (table < 0 || table >= tables()) {
563  throw GException::out_of_range(G_INX_OPERATOR2, "Table index",
564  table, tables());
565  }
566  #endif
567 
568  // Set indices and weighting factors for interpolation
569  update(arg1, arg2);
570 
571  // Perform 2D interpolation
572  double result = m_wgt1 * m_tables[table][m_inx1] +
576 
577  // Return result
578  return result;
579 }
580 
581 
582 /***********************************************************************//**
583  * @brief Trilinear interpolation operator for 3D tables
584  *
585  * @param[in] table Table index [0,...,tables()-1].
586  * @param[in] arg1 Value for first axis.
587  * @param[in] arg2 Value for second axis.
588  * @param[in] arg3 Value for second axis.
589  * @return Trilinearly interpolated response table.
590  *
591  * @exception GException::invalid_value
592  * Invalid response table dimension.
593  * @exception GException::out_of_range
594  * Invalid table index.
595  *
596  * Evaluates one response table at a given value. This method applies to
597  * three-dimensional response tables. If the table has a different dimension
598  * an exception is thrown.
599  *
600  * The evaluation is performed by a trilinear interpolation of the table
601  * values. If the specified value lies outside the range covered by the
602  * table, the table is linearly extrapolated from using either the first
603  * or the last two vector elements.
604  *
605  * @todo Write down formula.
606  ***************************************************************************/
607 double GCTAResponseTable::operator()(const int& table,
608  const double& arg1,
609  const double& arg2,
610  const double& arg3) const
611 {
612  // Throw exception if table is not 3D
613  if (axes() != 3) {
614  std::string msg = "Invalid response table dimension "+
615  gammalib::str(axes())+" encountered. Response "
616  "table needs to be three-dimensional.";
618  }
619 
620  // Throw exception if index is not valid
621  #if defined(G_RANGE_CHECK)
622  if (table < 0 || table >= tables()) {
623  throw GException::out_of_range(G_INX_OPERATOR3, "Table index",
624  table, tables());
625  }
626  #endif
627 
628  // Set indices and weighting factors for interpolation
629  update(arg1, arg2, arg3);
630 
631  // Perform 3D interpolation
632  double result = m_wgt1 * m_tables[table][m_inx1] +
640 
641  // Return result
642  return result;
643 }
644 
645 
646 /*==========================================================================
647  = =
648  = Public methods =
649  = =
650  ==========================================================================*/
651 
652 /***********************************************************************//**
653  * @brief Clear response table
654  *
655  * Clears the response table.
656  ***************************************************************************/
658 {
659  // Free class members
660  free_members();
661 
662  // Initialise members
663  init_members();
664 
665  // Return
666  return;
667 }
668 
669 
670 /***********************************************************************//**
671  * @brief Clone response table
672  *
673  * @return Deep copy of response table instance.
674  *
675  * Returns a pointer to a deep copy of the response table.
676  ***************************************************************************/
678 {
679  return new GCTAResponseTable(*this);
680 }
681 
682 
683 /***********************************************************************//**
684  * @brief Check whether a table exists
685  *
686  * @param[in] name Column name
687  * @return True if table exists.
688  *
689  * Checks whether a table exists.
690  ***************************************************************************/
691 bool GCTAResponseTable::has_table(const std::string& name) const
692 {
693  // Initialise existence flag
694  bool exists = false;
695 
696  // Loop over all tables to find index
697  int table = 0;
698  for (; table < tables(); ++table) {
699  if (m_colname_table[table] == name) {
700  exists = true;
701  break;
702  }
703  }
704 
705  // Return existence flag
706  return exists;
707 }
708 
709 
710 /***********************************************************************//**
711  * @brief Check whether an axis exists
712  *
713  * @param[in] name Column name
714  * @return True if axis exists.
715  *
716  * Check whether an axis exists.
717  ***************************************************************************/
718 bool GCTAResponseTable::has_axis(const std::string& name) const
719 {
720  // Initialise existence flag
721  bool exists = false;
722 
723  // Build column names
724  std::string col_lo = name + "_LO";
725  std::string col_hi = name + "_HI";
726 
727  // Loop over all axes to find index
728  int axis = 0;
729  for (; axis < axes(); ++axis) {
730  if ((axis_lo_name(axis) == col_lo) &&
731  (axis_hi_name(axis) == col_hi)) {
732  exists = true;
733  break;
734  }
735  }
736 
737  // Return existence flag
738  return exists;
739 }
740 
741 
742 /***********************************************************************//**
743  * @brief Determine index of table
744  *
745  * @param[in] name Column name
746  * @return Table index.
747  *
748  * @exception GException::invalid_value
749  * Table not found.
750  *
751  * Determines the index of a table @p name. An exception is thrown if the
752  * table is not found.
753  ***************************************************************************/
754 int GCTAResponseTable::table(const std::string& name) const
755 {
756  // Loop over all tables to find index
757  int table = 0;
758  for (; table < tables(); ++table) {
759  if (m_colname_table[table] == name) {
760  break;
761  }
762  }
763 
764  // Throw an exception if table has not been found
765  if (table >= tables()) {
766  std::string msg = "Table \""+name+"\" not found. Please verify the "
767  "response table.";
769  }
770 
771  // Return table index
772  return table;
773 }
774 
775 
776 /***********************************************************************//**
777  * @brief Return table unit
778  *
779  * @param[in] table Table index [0,...,tables()-1].
780  * @return Unit of table.
781  *
782  * @exception GException::out_of_range
783  * Table index out of range.
784  *
785  * Returns the unit of the table.
786  ***************************************************************************/
787 const std::string& GCTAResponseTable::unit(const int& table) const
788 {
789  // Optionally check if the table index is valid
790  #if defined(G_RANGE_CHECK)
791  if (table < 0 || table >= tables()) {
792  throw GException::out_of_range(G_UNIT, "Table index", table, tables());
793  }
794  #endif
795 
796  // Return units
797  return (m_units_table[table]);
798 }
799 
800 
801 /***********************************************************************//**
802  * @brief Scale table
803  *
804  * @param[in] table Table index [0,...,tables()-1].
805  * @param[in] scale Scaling factor.
806  *
807  * @exception GException::out_of_range
808  * Axis index out of range.
809  *
810  * Multiplies all values in a table by the scaling factor @p scale.
811  ***************************************************************************/
812 void GCTAResponseTable::scale(const int& table, const double& scale)
813 {
814  // Optionally check if the index is valid
815  #if defined(G_RANGE_CHECK)
816  if (table < 0 || table >= tables()) {
817  throw GException::out_of_range(G_SCALE, "Table index", table,
818  tables());
819  }
820  #endif
821 
822  // Scale table values
823  for (int i = 0; i < m_nelements; ++i) {
824  m_tables[table][i] *= scale;
825  }
826 
827  // Return
828  return;
829 }
830 
831 
832 /***********************************************************************//**
833  * @brief Determine index of an axis
834  *
835  * @param[in] name Column name
836  * @return Axis index.
837  *
838  * @exception GException::invalid_value
839  * Axis definition columns not found.
840  *
841  * Determines the index of an axis @p name, where the lower and upper bin
842  * edges of the axis are stored in the columns @p name_LO and @p name_HI.
843  * An exception is thrown if the axis is not found.
844  ***************************************************************************/
845 int GCTAResponseTable::axis(const std::string& name) const
846 {
847  // Build column names
848  std::string col_lo = name + "_LO";
849  std::string col_hi = name + "_HI";
850 
851  // Loop over all axes to find index
852  int axis = 0;
853  for (; axis < axes(); ++axis) {
854  if ((axis_lo_name(axis) == col_lo) &&
855  (axis_hi_name(axis) == col_hi)) {
856  break;
857  }
858  }
859 
860  // Throw an exception if axis has not been found
861  if (axis >= axes()) {
862  std::string msg = "Axis definition columns for axis \""+name+"\" not "
863  "found. Please verify the axis names in the "
864  "response table.";
865  throw GException::invalid_value(G_AXIS, msg);
866  }
867 
868  // Return axis index
869  return axis;
870 }
871 
872 
873 /***********************************************************************//**
874  * @brief Return number bins in an axis
875  *
876  * @param[in] axis Axis index [0,...,axes()-1].
877  * @return Number of bins along the specified axis.
878  *
879  * @exception GException::out_of_range
880  * Axis index out of range.
881  *
882  * Returns the number of bins in the specified @p axis.
883  ***************************************************************************/
884 int GCTAResponseTable::axis_bins(const int& axis) const
885 {
886  // Optionally check if the index is valid
887  #if defined(G_RANGE_CHECK)
888  if (axis < 0 || axis >= axes()) {
889  throw GException::out_of_range(G_AXIS_BINS, "Axis index", axis, axes());
890  }
891  #endif
892 
893  // Return axis length
894  return (m_axis_lo[axis].size());
895 }
896 
897 
898 /***********************************************************************//**
899  * @brief Return lower bin boundary for bin in axis
900  *
901  * @param[in] axis Axis index [0,...,axes()-1].
902  * @param[in] bin Bin index [0,...,axis_bins(axis)-1].
903  * @return Lower bin boundary.
904  *
905  * @exception GException::out_of_range
906  * Axis or bin index out of range.
907  *
908  * Returns the lower boundary for a given @p bin of the specified @p axis.
909  ***************************************************************************/
910 const double& GCTAResponseTable::axis_lo(const int& axis, const int& bin) const
911 {
912  // Check if index is valid
913  #if defined(G_RANGE_CHECK)
914  if (axis < 0 || axis >= axes()) {
915  throw GException::out_of_range(G_AXIS_LO, "Axis index", axis, axes());
916  }
917  if (bin < 0 || bin >= m_axis_lo[axis].size()) {
918  throw GException::out_of_range(G_AXIS_LO, "Bin index", bin,
919  m_axis_lo[axis].size());
920  }
921  #endif
922 
923  // Return bin boundary
924  return (m_axis_lo[axis][bin]);
925 }
926 
927 
928 /***********************************************************************//**
929  * @brief Return upper bin boundary for bin in axis
930  *
931  * @param[in] axis Axis index [0,...,axes()-1].
932  * @param[in] bin Bin index [0,...,axis_bins(axis)-1].
933  * @return Upper bin boundary.
934  *
935  * @exception GException::out_of_range
936  * Axis or bin index out of range.
937  *
938  * Returns the upper boundary for a given @p bin of the specified @p axis.
939  ***************************************************************************/
940 const double& GCTAResponseTable::axis_hi(const int& axis, const int& bin) const
941 {
942  // Optionally check if the index is valid
943  #if defined(G_RANGE_CHECK)
944  if (axis < 0 || axis >= axes()) {
945  throw GException::out_of_range(G_AXIS_HI, "Axis index", axis, axes());
946  }
947  if (bin < 0 || bin >= m_axis_hi[axis].size()) {
948  throw GException::out_of_range(G_AXIS_HI, "Bin index", bin,
949  m_axis_hi[axis].size());
950  }
951  #endif
952 
953  // Return bin boundary
954  return (m_axis_hi[axis][bin]);
955 }
956 
957 
958 /***********************************************************************//**
959  * @brief Return lower bin boundary FITS table column name for axis
960  *
961  * @param[in] axis Axis index [0,...,axes()-1].
962  * @return FITS table column name of lower bin boundary.
963  *
964  * @exception GException::out_of_range
965  * Axis index out of range.
966  *
967  * Returns the FITS table column name of the lower boundary for the specified
968  * @p axis.
969  ***************************************************************************/
970 const std::string& GCTAResponseTable::axis_lo_name(const int& axis) const
971 {
972  // Optionally check if the index is valid
973  #if defined(G_RANGE_CHECK)
974  if (axis < 0 || axis >= axes()) {
975  throw GException::out_of_range(G_AXIS_LO_NAME, "Axis index", axis, axes());
976  }
977  #endif
978 
979  // Return units
980  return (m_colname_lo[axis]);
981 }
982 
983 
984 /***********************************************************************//**
985  * @brief Return upper bin boundary FITS table column name for axis
986  *
987  * @param[in] axis Axis index [0,...,axes()-1].
988  * @return FITS table column name of upper bin boundary.
989  *
990  * @exception GException::out_of_range
991  * Axis index out of range.
992  *
993  * Returns the FITS table column name of the upper boundary for the specified
994  * @p axis.
995  ***************************************************************************/
996 const std::string& GCTAResponseTable::axis_hi_name(const int& axis) const
997 {
998  // Optionally check if the index is valid
999  #if defined(G_RANGE_CHECK)
1000  if (axis < 0 || axis >= axes()) {
1001  throw GException::out_of_range(G_AXIS_HI_NAME, "Axis index", axis, axes());
1002  }
1003  #endif
1004 
1005  // Return units
1006  return (m_colname_hi[axis]);
1007 }
1008 
1009 
1010 /***********************************************************************//**
1011  * @brief Return lower bin boundary unit for axis
1012  *
1013  * @param[in] axis Axis index [0,...,axes()-1].
1014  * @return Unit of lower bin boundary.
1015  *
1016  * @exception GException::out_of_range
1017  * Axis index out of range.
1018  *
1019  * Returns the unit of the lower boundary for the specified @p axis.
1020  ***************************************************************************/
1021 const std::string& GCTAResponseTable::axis_lo_unit(const int& axis) const
1022 {
1023  // Optionally check if the index is valid
1024  #if defined(G_RANGE_CHECK)
1025  if (axis < 0 || axis >= axes()) {
1026  throw GException::out_of_range(G_AXIS_LO_UNIT, "Axis index", axis, axes());
1027  }
1028  #endif
1029 
1030  // Return units
1031  return (m_units_lo[axis]);
1032 }
1033 
1034 
1035 /***********************************************************************//**
1036  * @brief Return upper bin boundary unit for axis
1037  *
1038  * @param[in] axis Axis index [0,...,axes()-1].
1039  * @return Unit of upper bin boundary.
1040  *
1041  * @exception GException::out_of_range
1042  * Axis index out of range.
1043  *
1044  * Returns the unit of the upper boundary for the specified @p axis.
1045  ***************************************************************************/
1046 const std::string& GCTAResponseTable::axis_hi_unit(const int& axis) const
1047 {
1048  // Optionally check if the index is valid
1049  #if defined(G_RANGE_CHECK)
1050  if (axis < 0 || axis >= axes()) {
1051  throw GException::out_of_range(G_AXIS_HI_UNIT, "Axis index", axis, axes());
1052  }
1053  #endif
1054 
1055  // Return units
1056  return (m_units_hi[axis]);
1057 }
1058 
1059 
1060 /***********************************************************************//**
1061  * @brief Set nodes for a linear axis
1062  *
1063  * @param[in] axis Axis index [0,...,axes()-1].
1064  *
1065  * @exception GException::out_of_range
1066  * Axis index out of range.
1067  *
1068  * Set axis nodes so that each node is the linear mean of the lower and upper
1069  * bin boundary, i.e.
1070  * \f[ n_i = 0.5 \times ({\rm LO}_i + {\rm HI}_i) \f]
1071  * where
1072  * \f$n_i\f$ is node \f$i\f$,
1073  * \f${\rm LO}_i\f$ is the lower bin boundary for bin \f$i\f$, and
1074  * \f${\rm HI}_i\f$ is the upper bin boundary for bin \f$i\f$.
1075  ***************************************************************************/
1076 void GCTAResponseTable::axis_linear(const int& axis)
1077 {
1078  // Optionally check if the index is valid
1079  #if defined(G_RANGE_CHECK)
1080  if (axis < 0 || axis >= axes()) {
1081  throw GException::out_of_range(G_AXIS_LINEAR, "Axis index", axis,
1082  axes());
1083  }
1084  #endif
1085 
1086  // Get number of bins
1087  int bins = axis_bins(axis);
1088 
1089  // Allocate node values
1090  std::vector<double> axis_nodes(bins);
1091 
1092  // Compute nodes
1093  for (int i = 0; i < bins; ++i) {
1094  axis_nodes[i] = 0.5*(m_axis_lo[axis][i] + m_axis_hi[axis][i]);
1095  }
1096 
1097  // Set node array
1098  m_axis_nodes[axis] = GNodeArray(axis_nodes);
1099 
1100  // Return
1101  return;
1102 }
1103 
1104 
1105 /***********************************************************************//**
1106  * @brief Set nodes for a logarithmic (base 10) axis
1107  *
1108  * @param[in] axis Axis index [0,...,axes()-1].
1109  *
1110  * @exception GException::out_of_range
1111  * Axis index out of range.
1112  *
1113  * Set axis nodes so that each node is the logarithmic mean of the lower and
1114  * upper bin boundary, i.e.
1115  * \f[ n_i = \log \sqrt{{\rm LO}_i \times {\rm HI}_i} \f]
1116  * where
1117  * \f$n_i\f$ is node \f$i\f$,
1118  * \f${\rm LO}_i\f$ is the lower bin boundary for bin \f$i\f$, and
1119  * \f${\rm HI}_i\f$ is the upper bin boundary for bin \f$i\f$.
1120  *
1121  * @todo Check that none of the axis boundaries is non-positive.
1122  ***************************************************************************/
1123 void GCTAResponseTable::axis_log10(const int& axis)
1124 {
1125  // Optionally check if the index is valid
1126  #if defined(G_RANGE_CHECK)
1127  if (axis < 0 || axis >= axes()) {
1128  throw GException::out_of_range(G_AXIS_LOG10, "Axis index", axis,
1129  axes());
1130  }
1131  #endif
1132 
1133  // Get number of bins
1134  int bins = axis_bins(axis);
1135 
1136  // Allocate node values
1137  std::vector<double> axis_nodes(bins);
1138 
1139  // Compute nodes
1140  for (int i = 0; i < bins; ++i) {
1141  axis_nodes[i] = std::log10(std::sqrt(m_axis_lo[axis][i] *
1142  m_axis_hi[axis][i]));
1143  }
1144 
1145  // Set node array
1146  m_axis_nodes[axis] = GNodeArray(axis_nodes);
1147 
1148  // Return
1149  return;
1150 }
1151 
1152 
1153 /***********************************************************************//**
1154  * @brief Set nodes for a radians axis
1155  *
1156  * @param[in] axis Axis index [0,...,axes()-1].
1157  *
1158  * @exception GException::out_of_range
1159  * Axis index out of range.
1160  *
1161  * Set axis nodes so that each node is the lo mean of the lower and upper
1162  * bin boundary in radians, i.e.
1163  * \f[ n_i = \frac{\pi}{360} \times ({\rm LO}_i + {\rm HI}_i) \f]
1164  * where
1165  * \f$n_i\f$ is node \f$i\f$,
1166  * \f${\rm LO}_i\f$ is the lower bin boundary for bin \f$i\f$, and
1167  * \f${\rm HI}_i\f$ is the upper bin boundary for bin \f$i\f$.
1168  ***************************************************************************/
1169 void GCTAResponseTable::axis_radians(const int& axis)
1170 {
1171  // Optionally check if the index is valid
1172  #if defined(G_RANGE_CHECK)
1173  if (axis < 0 || axis >= axes()) {
1174  throw GException::out_of_range(G_AXIS_RADIANS, "Axis index", axis,
1175  axes());
1176  }
1177  #endif
1178 
1179  // Get number of bins
1180  int bins = axis_bins(axis);
1181 
1182  // Allocate node values
1183  std::vector<double> axis_nodes(bins);
1184 
1185  // Compute nodes
1186  for (int i = 0; i < bins; ++i) {
1187  axis_nodes[i] = 0.5*(m_axis_lo[axis][i] + m_axis_hi[axis][i]) *
1189  }
1190 
1191  // Set node array
1192  m_axis_nodes[axis] = GNodeArray(axis_nodes);
1193 
1194  // Return
1195  return;
1196 }
1197 
1198 
1199 /***********************************************************************//**
1200  * @brief Append an axis to the response table
1201  *
1202  * @param[in] axis_lo Lower axis boundaries.
1203  * @param[in] axis_hi Upper axis boundaries.
1204  * @param[in] name Axis name.
1205  * @param[in] unit Axis unit.
1206  *
1207  * @exception GException::invalid_argument
1208  * Incompatible axis boundaries.
1209  *
1210  * Append an axis to the response table.
1211  ***************************************************************************/
1212 void GCTAResponseTable::append_axis(const std::vector<double>& axis_lo,
1213  const std::vector<double>& axis_hi,
1214  const std::string& name,
1215  const std::string& unit)
1216 {
1217  // Throw an exception if the length of axis_lo and axis_hi are different
1218  if (axis_lo.size() != axis_hi.size()) {
1219  std::string msg = "Number of elements in lower axis boundaries "
1220  "mismatches the number of elements in upper axis "
1221  "boundaries ("+gammalib::str(axis_lo.size())+
1222  " != "+gammalib::str(axis_hi.size())+").";
1224  }
1225 
1226  // Set axis names
1227  std::string name_lo = name + "_LO";
1228  std::string name_hi = name + "_HI";
1229 
1230  // Append axis
1231  m_colname_lo.push_back(name_lo);
1232  m_colname_hi.push_back(name_hi);
1233  m_axis_lo.push_back(axis_lo);
1234  m_axis_hi.push_back(axis_hi);
1235  m_units_lo.push_back(unit);
1236  m_units_hi.push_back(unit);
1237 
1238 
1239  // Set node array
1240  std::vector<double> axis_nodes(axis_lo.size());
1241  for (int k = 0; k < axis_lo.size(); ++k) {
1242  axis_nodes[k] = 0.5*(axis_lo[k] + axis_hi[k]);
1243  }
1244  m_axis_nodes.push_back(GNodeArray(axis_nodes));
1245 
1246  // Increment number of axes
1247  m_naxes++;
1248 
1249  // Compute the cube size
1250  m_nelements = axis_bins(0);
1251  for (int i = 1; i < axes(); ++i) {
1252  m_nelements *= axis_bins(i);
1253  }
1254 
1255  // Return
1256  return;
1257 }
1258 
1259 
1260 /***********************************************************************//**
1261  * @brief Append table to response table
1262  *
1263  * @param[in] name Table name.
1264  * @param[in] unit Table unit.
1265  *
1266  * @exception GException::invalid_value
1267  * No axes have been specified.
1268  *
1269  * Append a table to the response table. The number of elements in the table
1270  * is the product of the length of all axes. All elements are set to zero by
1271  * default.
1272  ***************************************************************************/
1273 void GCTAResponseTable::append_table(const std::string& name,
1274  const std::string& unit)
1275 {
1276  // Throw an exception message when m_nelement is zero
1277  if (m_nelements == 0) {
1278  std::string msg = "No axis columns have been defined. Please append "
1279  "axis columns before appending table columns.";
1281  }
1282 
1283  // Append table column name and unit
1284  m_colname_table.push_back(name);
1285  m_units_table.push_back(unit);
1286 
1287  // Initialise empty table column
1288  std::vector<double> table(m_nelements, 0.0);
1289 
1290  // Append column
1291  m_tables.push_back(table);
1292 
1293  // Increment number of table columns
1294  m_ntables++;
1295 
1296  // Return
1297  return;
1298 }
1299 
1300 
1301 /***********************************************************************//**
1302  * @brief Return axis nodes
1303  *
1304  * @param[in] axis Axis index [0,...,axes()-1].
1305  * @return Node array for axis.
1306  *
1307  * @exception GException::out_of_range
1308  * Axis index out of range.
1309  *
1310  * Returns the node array of the specified @p axis.
1311  ***************************************************************************/
1312 const GNodeArray& GCTAResponseTable::axis_nodes(const int& axis) const
1313 {
1314  // Optionally check if the index is valid
1315  #if defined(G_RANGE_CHECK)
1316  if (axis < 0 || axis >= axes()) {
1317  throw GException::out_of_range(G_AXIS_NODES, "Axis index", axis,
1318  axes());
1319  }
1320  #endif
1321 
1322  // Return node array
1323  return (m_axis_nodes[axis]);
1324 }
1325 
1326 
1327 /***********************************************************************//**
1328  * @brief Read response table from FITS table HDU
1329  *
1330  * @param[in] table FITS table.
1331  *
1332  * Reads CTA response table information from a FITS table. The FITS table
1333  * is expected to have a single row, and axes and table information are
1334  * found in vector columns. Axes information is expected to be placed
1335  * before table information.
1336  *
1337  * Each axis is defined by two vector columns of equal width, describing the
1338  * lower and upper limits for each bin. The column names for this boundary
1339  * information terminates by "_LO" and "HI", respectively (upper case). It
1340  * is expected that the "_LO" column preceeds the "_HI" column.
1341  *
1342  * Following the axes columns are table columns of equal width. The width
1343  * of each table column is given by the product of the lengths of all
1344  * axes. It is furthermore expected that the first axis is the most rapidely
1345  * varying index of the vector.
1346  ***************************************************************************/
1348 {
1349  // Clear instance
1350  clear();
1351 
1352  // Read column names
1353  read_colnames(table);
1354 
1355  // Read axes
1356  read_axes(table);
1357 
1358  // Read tables
1359  read_tables(table);
1360 
1361  // Read TELESCOP keyword
1362  m_telescope = (table.has_card("TELESCOP")) ? table.string("TELESCOP") : "";
1363 
1364  // Return
1365  return;
1366 }
1367 
1368 
1369 /***********************************************************************//**
1370  * @brief Write response table into FITS table HDU
1371  *
1372  * @param[in] table FITS table.
1373  *
1374  * Writes the response table in a FITS table HDU.
1375  *
1376  * Axes columns are written in floating point columns with the suffixes "_LO"
1377  * and "_HI". After the axes columns the response tables are written in
1378  * floating point columns.
1379  ***************************************************************************/
1381 {
1382  // Initialise dimension vector
1383  std::vector<int> dim;
1384 
1385  // Loop over all response table axes
1386  for (int iaxis = 0; iaxis < m_naxes; ++iaxis) {
1387 
1388  // Create axis columns
1389  GFitsTableFloatCol col_lo(m_colname_lo[iaxis], 1,
1390  m_axis_lo[iaxis].size());
1391  GFitsTableFloatCol col_hi(m_colname_hi[iaxis], 1,
1392  m_axis_hi[iaxis].size());
1393 
1394  // Loop through all elements in this axis column
1395  for (int i = 0; i < m_axis_lo[iaxis].size(); ++i) {
1396  col_lo(0,i) = m_axis_lo[iaxis][i];
1397  col_hi(0,i) = m_axis_hi[iaxis][i];
1398  }
1399 
1400  // Set column units
1401  col_lo.unit(m_units_lo[iaxis]);
1402  col_hi.unit(m_units_hi[iaxis]);
1403 
1404  // Append column to FITS table
1405  table.append(col_lo);
1406  table.append(col_hi);
1407 
1408  // Append dimension
1409  dim.push_back(axis_bins(iaxis));
1410 
1411  } // endif: looped over all axes in response table
1412 
1413  // Loop over all tables
1414  for (int itable = 0; itable < m_ntables; ++itable) {
1415 
1416  // Create table column
1417  GFitsTableFloatCol col_table(m_colname_table[itable], 1,
1418  m_tables[itable].size());
1419 
1420  // Loop through elements in this table column
1421  for (int i = 0; i < m_tables[itable].size() ; ++i) {
1422  col_table(0,i) = m_tables[itable][i];
1423  }
1424 
1425  // Set column unit
1426  col_table.unit(m_units_table[itable]);
1427 
1428  // Set column dimension
1429  col_table.dim(dim);
1430 
1431  // Append column to table
1432  table.append(col_table);
1433 
1434  } // endfor: looped over all axes in response table
1435 
1436  // Write TELESCOPE keyword
1437  table.card("TELESCOP", m_telescope, "Telescope");
1438 
1439  // Return
1440  return;
1441 }
1442 
1443 
1444 /***********************************************************************//**
1445  * @brief Print response table information
1446  *
1447  * @param[in] chatter Chattiness.
1448  * @return String containing response table information.
1449  *
1450  * Puts CTA response table information into a std::string object for
1451  * printing.
1452  ***************************************************************************/
1453 std::string GCTAResponseTable::print(const GChatter& chatter) const
1454 {
1455  // Initialise result string
1456  std::string result;
1457 
1458  // Continue only if chatter is not silent
1459  if (chatter != SILENT) {
1460 
1461  // Append header
1462  result.append("=== GCTAResponseTable ===");
1463 
1464  // Append information
1465  result.append("\n"+gammalib::parformat("Telescope") +
1466  m_telescope);
1467  result.append("\n"+gammalib::parformat("Dimension") +
1468  gammalib::str(axes()));
1469 
1470  // Append axes information
1471  for (int i = 0; i < axes(); ++i) {
1472  result.append("\n"+gammalib::parformat("Axis " +
1473  gammalib::str(i)));
1474  result.append(gammalib::str(axis_bins(i)));
1475  result.append(" (");
1476  result.append(m_colname_lo[i]);
1477  if (!m_units_lo[i].empty()) {
1478  result.append(" ["+m_units_lo[i]+"]");
1479  }
1480  result.append(", ");
1481  result.append(m_colname_hi[i]);
1482  if (!m_units_hi[i].empty()) {
1483  result.append(" ["+m_units_hi[i]+"]");
1484  }
1485  result.append(")");
1486  }
1487 
1488  // Append table information
1489  for (int i = 0; i < tables(); ++i) {
1490  result.append("\n"+gammalib::parformat("Table " +
1491  gammalib::str(i)));
1492  result.append(m_colname_table[i]);
1493  if (!m_units_table[i].empty()) {
1494  result.append(" ["+m_units_table[i]+"]");
1495  }
1496  }
1497 
1498  } // endif: chatter was not silent
1499 
1500  // Return result
1501  return result;
1502 }
1503 
1504 
1505 /*==========================================================================
1506  = =
1507  = Private methods =
1508  = =
1509  ==========================================================================*/
1510 
1511 /***********************************************************************//**
1512  * @brief Initialise response table members
1513  *
1514  * Initialises all members of the response table.
1515  ***************************************************************************/
1517 {
1518  // Initialise members
1519  m_naxes = 0;
1520  m_ntables = 0;
1521  m_nelements = 0;
1522  m_colname_lo.clear();
1523  m_colname_hi.clear();
1524  m_colname_table.clear();
1525  m_axis_lo.clear();
1526  m_axis_hi.clear();
1527  m_units_lo.clear();
1528  m_units_hi.clear();
1529  m_units_table.clear();
1530  m_axis_nodes.clear();
1531  m_tables.clear();
1532  m_telescope.clear();
1533 
1534  // Initialise cache
1535  m_inx_left = 0;
1536  m_inx_right = 0;
1537  m_wgt_left = 0.0;
1538  m_wgt_right = 0.0;
1539  m_inx1 = 0;
1540  m_inx2 = 0;
1541  m_inx3 = 0;
1542  m_inx4 = 0;
1543  m_inx5 = 0;
1544  m_inx6 = 0;
1545  m_inx7 = 0;
1546  m_inx8 = 0;
1547  m_wgt1 = 0.0;
1548  m_wgt2 = 0.0;
1549  m_wgt3 = 0.0;
1550  m_wgt4 = 0.0;
1551  m_wgt5 = 0.0;
1552  m_wgt6 = 0.0;
1553  m_wgt7 = 0.0;
1554  m_wgt8 = 0.0;
1555 
1556  // Return
1557  return;
1558 }
1559 
1560 
1561 /***********************************************************************//**
1562  * @brief Copy response table members
1563  *
1564  * @param[in] table Response table.
1565  *
1566  * Copies all response table members.
1567  ***************************************************************************/
1569 {
1570  // Copy number of bins
1571  m_naxes = table.m_naxes;
1572  m_ntables = table.m_ntables;
1573  m_nelements = table.m_nelements;
1574  m_colname_lo = table.m_colname_lo;
1575  m_colname_hi = table.m_colname_hi;
1577  m_axis_lo = table.m_axis_lo;
1578  m_axis_hi = table.m_axis_hi;
1579  m_units_lo = table.m_units_lo;
1580  m_units_hi = table.m_units_hi;
1581  m_units_table = table.m_units_table;
1582  m_axis_nodes = table.m_axis_nodes;
1583  m_tables = table.m_tables;
1584  m_telescope = table.m_telescope;
1585 
1586  // Copy cache
1587  m_inx_left = table.m_inx_left;
1588  m_inx_right = table.m_inx_right;
1589  m_wgt_left = table.m_wgt_left;
1590  m_wgt_right = table.m_wgt_right;
1591  m_inx1 = table.m_inx1;
1592  m_inx2 = table.m_inx2;
1593  m_inx3 = table.m_inx3;
1594  m_inx4 = table.m_inx4;
1595  m_inx5 = table.m_inx5;
1596  m_inx6 = table.m_inx6;
1597  m_inx7 = table.m_inx7;
1598  m_inx8 = table.m_inx8;
1599  m_wgt1 = table.m_wgt1;
1600  m_wgt2 = table.m_wgt2;
1601  m_wgt3 = table.m_wgt3;
1602  m_wgt4 = table.m_wgt4;
1603  m_wgt5 = table.m_wgt5;
1604  m_wgt6 = table.m_wgt6;
1605  m_wgt7 = table.m_wgt7;
1606  m_wgt8 = table.m_wgt8;
1607 
1608  // Return
1609  return;
1610 }
1611 
1612 
1613 /***********************************************************************//**
1614  * @brief Delete response table members
1615  *
1616  * De-allocates any memory that was allocated by the response table.
1617  ***************************************************************************/
1619 {
1620  // Return
1621  return;
1622 }
1623 
1624 
1625 /***********************************************************************//**
1626  * @brief Read column names from FITS HDU
1627  *
1628  * @param[in] hdu Response table HDU.
1629  *
1630  * @exception GException::invalid_value
1631  * Invalid response table format encountered.
1632  *
1633  * Read the response table column names from the HDU. Column names
1634  * terminating with "_LO" and "_HI" define the table axes, while all other
1635  * column names define response tables. It is assumed that table axes are
1636  * given in subsequent order, with the first column corresponding to the
1637  * lower bin boundaries of the first axis. Lower bin boundaries are indicated
1638  * by the "_LO" termination. Upper bin boundaries are assumed to follow
1639  * immediately the lower bin boundaires and are designated by the "_HI"
1640  * termination.
1641  *
1642  * This method sets the following members:
1643  * m_colname_lo - Column names of lower boundaries
1644  * m_colname_hi - Column names of upper boundaries
1645  * m_colname_table - Column names of tables
1646  * m_naxes - Number of axes
1647  * m_ntables - Number of tables
1648  *
1649  * @todo Implement exceptions for invalid HDU format
1650  ***************************************************************************/
1652 {
1653  // Clear column name arrays
1654  m_naxes = 0;
1655  m_ntables = 0;
1656  m_colname_lo.clear();
1657  m_colname_hi.clear();
1658  m_colname_table.clear();
1659 
1660  // Initialise search mode. There are three search modes:
1661  // 0 - we're looking for the next axis by searching for a column
1662  // terminating with "_LO"
1663  // 1 - we're looking for the upper boundary of an axis, terminating
1664  // with "_HI"
1665  // 2 - we're looking for a table column
1666  int mode = 0;
1667  std::string lo_column;
1668  std::string next_column;
1669 
1670  // Extract column names for all axes
1671  for (int i = 0; i < hdu.ncols(); ++i) {
1672 
1673  // Get column name and unit
1674  std::string colname = hdu[i]->name();
1675 
1676  // If we search for a "_LO" column, check if we have one. If one
1677  // is found, change the search mode to 1 and set the expected name
1678  // for the "_HI" column. If none is found, change the search
1679  // mode to 2 since from now on we should only have table
1680  // columns.
1681  if (mode == 0) {
1682  size_t pos = colname.rfind("_LO");
1683  if (pos != std::string::npos) {
1684  mode = 1;
1685  lo_column = colname;
1686  next_column = colname.substr(0, pos) + "_HI";
1687  }
1688  else {
1689  if (colname.rfind("_HI") != std::string::npos) {
1690  std::string msg = "Column \""+colname+"\" encountered "
1691  "without preceeding \"_LO\" column. "
1692  "Please correct response table format.";
1694  }
1695  else {
1696  mode = 2;
1697  m_colname_table.push_back(colname);
1698  }
1699  }
1700  }
1701 
1702  // If we search for a "_HI" column, check if we have the
1703  // expected column name. If this is the case, switch back to
1704  // search mode 0 to get the next "_LO" column. Otherwise
1705  // throw an exception.
1706  else if (mode == 1) {
1707  if (colname == next_column) {
1708  mode = 0;
1709  m_colname_lo.push_back(lo_column);
1710  m_colname_hi.push_back(next_column);
1711  }
1712  else {
1713  std::string msg = "Expected column \""+next_column+"\" not "
1714  "found. The axis upper bin boundary "
1715  "column has to follow immediately the "
1716  "lower bin boundary column.";
1718  }
1719  }
1720 
1721  // If we search for a table column, make sure that we have
1722  // neither a "_LO" nor a "_HI" column.
1723  else {
1724  if (colname.rfind("_LO") != std::string::npos) {
1725  std::string msg = "Column \""+colname+"\" encountered while "
1726  "searching for table columns. All "
1727  "lower bin boundary columns have to be "
1728  "placed before the table columns.";
1730  }
1731  else if (colname.rfind("_HI") != std::string::npos) {
1732  std::string msg = "Column \""+colname+"\" encountered while "
1733  "searching for table columns. All "
1734  "upper bin boundary columns have to be "
1735  "placed before the table columns.";
1737  }
1738  else {
1739  m_colname_table.push_back(colname);
1740  }
1741  }
1742 
1743  } // endfor: looped over all columns
1744 
1745  // Store number of axes
1746  m_naxes = m_colname_lo.size();
1747 
1748  // Store number of tables
1749  m_ntables = m_colname_table.size();
1750 
1751  // Return
1752  return;
1753 }
1754 
1755 
1756 /***********************************************************************//**
1757  * @brief Read axes definitions from FITS HDU
1758  *
1759  * @param[in] hdu Response table HDU.
1760  *
1761  * @exception GException::invalid_value
1762  * Incompatible axes columns encoutered.
1763  *
1764  * Reads the lower and upper boundaries for all response table axes from the
1765  * FITS HDU. The method verifies if the lower and upper boundary axes have
1766  * the same vector dimensions.
1767  *
1768  * The method also allocates the nodes for each axis, assuming that axis will
1769  * be used linearily, where each node is given by the linear mean of the
1770  * lower and upper boundary. In case that the axis should be used
1771  * logarithmically, the method axis_log10() can be used to change to a
1772  * logarithmic scale. The axis_linear() can be used to switch back to a
1773  * linear scale.
1774  *
1775  * This method sets the following members:
1776  * m_axis_lo - Axes lower boundaries
1777  * m_axis_hi - Axes upper boundaries
1778  * m_axis_nodes - Axes mean values
1779  ***************************************************************************/
1781 {
1782  // Clear axes arrays
1783  m_axis_lo.clear();
1784  m_axis_hi.clear();
1785  m_axis_nodes.clear();
1786 
1787  // Loop over all dimensions
1788  for (int i = 0; i < axes(); ++i) {
1789 
1790  // Get pointers to table columns
1791  const GFitsTableCol* col_lo = hdu[m_colname_lo[i]];
1792  const GFitsTableCol* col_hi = hdu[m_colname_hi[i]];
1793 
1794  // Extract number of bins. Make sure that both columns have the
1795  // same number of bins
1796  int num = col_lo->number();
1797  if (num != col_hi->number()) {
1798  std::string msg = "Axis lower bin boundary column \""+
1799  m_colname_lo[i]+"\" contains "+
1800  gammalib::str(num)+" elements while upper "
1801  "bin boundary column \""+m_colname_hi[i]+
1802  "\" contains "+gammalib::str(col_hi->number())+
1803  " elements. Both columns need to contain the "
1804  "same number of elements.";
1806  }
1807 
1808  // Initialise axis and node arrays
1809  std::vector<double> axis_lo(num);
1810  std::vector<double> axis_hi(num);
1811  std::vector<double> axis_nodes(num);
1812 
1813  // Copy axis information into arrays
1814  for (int k = 0; k < num; ++k) {
1815  axis_lo[k] = col_lo->real(0,k);
1816  axis_hi[k] = col_hi->real(0,k);
1817  axis_nodes[k] = 0.5*(axis_lo[k] + axis_hi[k]);
1818  }
1819 
1820  // Push axis array on storage
1821  m_axis_lo.push_back(axis_lo);
1822  m_axis_hi.push_back(axis_hi);
1823 
1824  // Push units on storage
1825  m_units_lo.push_back(col_lo->unit());
1826  m_units_hi.push_back(col_hi->unit());
1827 
1828  // Create node array
1829  m_axis_nodes.push_back(GNodeArray(axis_nodes));
1830 
1831  } // endfor: looped over all dimensions
1832 
1833  // Return
1834  return;
1835 }
1836 
1837 
1838 /***********************************************************************//**
1839  * @brief Read tables
1840  *
1841  * @param[in] hdu Response table HDU.
1842  *
1843  * @exception GException::invalid_value
1844  * Table vector of bad size encountered.
1845  *
1846  * Reads the tables from the FITS HDU.
1847  *
1848  * The method also sets m_nelements which gives the number of elements in a
1849  * response table.
1850  ***************************************************************************/
1852 {
1853  // Clear tables
1854  m_tables.clear();
1855 
1856  // Compute expected cube size
1857  m_nelements = (axes() > 0) ? axis_bins(0) : 0;
1858  for (int i = 1; i < axes(); ++i) {
1859  m_nelements *= axis_bins(i);
1860  }
1861 
1862  // Loop over all tables
1863  for (int i = 0; i < tables(); ++i) {
1864 
1865  // Get pointer to table column
1866  const GFitsTableCol* col = hdu[m_colname_table[i]];
1867 
1868  // Extract number of bins. Verify that the number of bins
1869  // corresponds to the expectation.
1870  int num = col->number();
1871  if (num != m_nelements) {
1872  std::string msg = "Table column \""+m_colname_table[i]+"\" "
1873  "contains "+gammalib::str(num)+" elements "
1874  "while "+gammalib::str(m_nelements)+" "
1875  "elements were expected from the axis "
1876  "definitions. Please check the response "
1877  "table for consistency.";
1879  }
1880 
1881  // Check table dimension if there is dimension information
1882  if (col->dim().size() > 0) {
1883 
1884  // Check that the table dimension is identical to the number of
1885  // axis pairs
1886  if (axes() != col->dim().size()) {
1887  std::string msg = "Table column \""+m_colname_table[i]+"\" "
1888  "dimension "+
1889  gammalib::str(col->dim().size())+" is "
1890  "inconsistent with "+
1891  gammalib::str(axes())+" axis pairs. Please "
1892  "check the response table for consistency.";
1894  }
1895 
1896  // Check for all axes that the length of the axis columns is
1897  // identical to the size of the table in that axis direction
1898  for (int i = 0; i < axes(); ++i) {
1899  if (axis_bins(i) != col->dim()[i]) {
1900  std::string msg = "Table column \""+m_colname_table[i]+"\" "
1901  "has "+gammalib::str(col->dim()[i])+
1902  " bins in axis "+gammalib::str(i)+
1903  " while corresponding axis columns \""+
1904  m_colname_lo[i]+"\" and \""+
1905  m_colname_hi[i]+"\" have a length of "+
1906  gammalib::str(axis_bins(i))+". Please "
1907  "check the response table for "
1908  "consistency.";
1910  }
1911  }
1912 
1913  } // endif: table had dimension information
1914 
1915 
1916  // Initialise table
1917  std::vector<double> table(num);
1918 
1919  // Copy table values
1920  for (int k = 0; k < num; ++k) {
1921  table[k] = col->real(0,k);
1922  }
1923 
1924  // Push table into storage
1925  m_tables.push_back(table);
1926 
1927  // Push units on storage
1928  m_units_table.push_back(col->unit());
1929 
1930  } // endfor: looped over all tables
1931 
1932  // Return
1933  return;
1934 }
1935 
1936 
1937 /***********************************************************************//**
1938  * @brief Update 1D cache
1939  *
1940  * @param[in] arg Argument.
1941  *
1942  * Updates the 1D interpolation cache. The interpolation cache is composed
1943  * of two indices and weights that define 2 data values of the 2D table
1944  * that are used for linear interpolation.
1945  *
1946  * @todo Write down formula
1947  ***************************************************************************/
1948 void GCTAResponseTable::update(const double& arg) const
1949 {
1950  // Get pointer to node array
1951  const GNodeArray* nodes = &(m_axis_nodes[0]);
1952 
1953  // Set value for node array
1954  nodes->set_value(arg);
1955 
1956  // Set indices and weighting factors for interpolation
1957  m_inx_left = nodes->inx_left();
1958  m_inx_right = nodes->inx_right();
1959  m_wgt_left = nodes->wgt_left();
1960  m_wgt_right = nodes->wgt_right();
1961 
1962  // Return
1963  return;
1964 }
1965 
1966 
1967 /***********************************************************************//**
1968  * @brief Update 2D cache
1969  *
1970  * @param[in] arg1 Argument for first axis.
1971  * @param[in] arg2 Argument for second axis.
1972  *
1973  * Updates the 2D interpolation cache. The interpolation cache is composed
1974  * of four indices and weights that define 4 data values of the 2D table
1975  * that are used for bilinear interpolation.
1976  *
1977  * @todo Write down formula
1978  ***************************************************************************/
1979 void GCTAResponseTable::update(const double& arg1, const double& arg2) const
1980 {
1981  // Get pointers to node arrays
1982  const GNodeArray* nodes1 = &(m_axis_nodes[0]);
1983  const GNodeArray* nodes2 = &(m_axis_nodes[1]);
1984 
1985  // Set values for node arrays
1986  nodes1->set_value(arg1);
1987  nodes2->set_value(arg2);
1988 
1989  // Compute offsets
1990  int size1 = axis_bins(0);
1991  int offset_left = nodes2->inx_left() * size1;
1992  int offset_right = nodes2->inx_right() * size1;
1993 
1994  // Set indices for bi-linear interpolation
1995  m_inx1 = nodes1->inx_left() + offset_left;
1996  m_inx2 = nodes1->inx_left() + offset_right;
1997  m_inx3 = nodes1->inx_right() + offset_left;
1998  m_inx4 = nodes1->inx_right() + offset_right;
1999 
2000  // Set weighting factors for bi-linear interpolation
2001  m_wgt1 = nodes1->wgt_left() * nodes2->wgt_left();
2002  m_wgt2 = nodes1->wgt_left() * nodes2->wgt_right();
2003  m_wgt3 = nodes1->wgt_right() * nodes2->wgt_left();
2004  m_wgt4 = nodes1->wgt_right() * nodes2->wgt_right();
2005 
2006  // Return
2007  return;
2008 }
2009 
2010 
2011 /***********************************************************************//**
2012  * @brief Update 3D cache
2013  *
2014  * @param[in] arg1 Argument for first axis.
2015  * @param[in] arg2 Argument for second axis.
2016  * @param[in] arg3 Argument for third axis.
2017  *
2018  * Updates the 3D interpolation cache. The interpolation cache is composed
2019  * of four indices and weights that define 4 data values of the 2D table
2020  * that are used for bilinear interpolation.
2021  *
2022  * @todo Write down formula
2023  ***************************************************************************/
2024 void GCTAResponseTable::update(const double& arg1, const double& arg2,
2025  const double& arg3) const
2026 {
2027 
2028  // Get pointers to node arrays (circumvent const correctness)
2029  const GNodeArray* nodes1 = &(m_axis_nodes[0]);
2030  const GNodeArray* nodes2 = &(m_axis_nodes[1]);
2031  const GNodeArray* nodes3 = &(m_axis_nodes[2]);
2032 
2033  // Set values for node arrays
2034  nodes1->set_value(arg1);
2035  nodes2->set_value(arg2);
2036  nodes3->set_value(arg3);
2037 
2038  // Compute offsets
2039  int size1 = axis_bins(0);
2040  int size2 = axis_bins(1);
2041  int size12 = size1 * size2;
2042  int offset_left_2 = nodes2->inx_left() * size1;
2043  int offset_right_2 = nodes2->inx_right() * size1;
2044  int offset_left_3 = nodes3->inx_left() * size12;
2045  int offset_right_3 = nodes3->inx_right() * size12;
2046 
2047  // Pre-compute stuff
2048  int inx1l2l = nodes1->inx_left() + offset_left_2;
2049  int inx1l2r = nodes1->inx_left() + offset_right_2;
2050  int inx1r2l = nodes1->inx_right() + offset_left_2;
2051  int inx1r2r = nodes1->inx_right() + offset_right_2;
2052 
2053  // Set indices for tri-linear interpolation
2054  m_inx1 = inx1l2l + offset_left_3;
2055  m_inx2 = inx1l2l + offset_right_3;
2056  m_inx3 = inx1l2r + offset_left_3;
2057  m_inx4 = inx1l2r + offset_right_3;
2058  m_inx5 = inx1r2l + offset_left_3;
2059  m_inx6 = inx1r2l + offset_right_3;
2060  m_inx7 = inx1r2r + offset_left_3;
2061  m_inx8 = inx1r2r + offset_right_3;
2062 
2063  // Pre-compute stuff
2064  double wgt1l2l = nodes1->wgt_left() * nodes2->wgt_left();
2065  double wgt1l2r = nodes1->wgt_left() * nodes2->wgt_right();
2066  double wgt1r2l = nodes1->wgt_right() * nodes2->wgt_left();
2067  double wgt1r2r = nodes1->wgt_right() * nodes2->wgt_right();
2068 
2069  // Set weighting factors for tri-linear interpolation
2070  m_wgt1 = wgt1l2l * nodes3->wgt_left();
2071  m_wgt2 = wgt1l2l * nodes3->wgt_right();
2072  m_wgt3 = wgt1l2r * nodes3->wgt_left();
2073  m_wgt4 = wgt1l2r * nodes3->wgt_right();
2074  m_wgt5 = wgt1r2l * nodes3->wgt_left();
2075  m_wgt6 = wgt1r2l * nodes3->wgt_right();
2076  m_wgt7 = wgt1r2r * nodes3->wgt_left();
2077  m_wgt8 = wgt1r2r * nodes3->wgt_right();
2078 
2079  // Return
2080  return;
2081 }
int m_inx_left
Index of left node.
const GNodeArray & axis_nodes(const int &axis) const
Return axis nodes.
void unit(const std::string &unit)
Set column unit.
const int & ncols(void) const
Return number of columns in table.
Definition: GFitsTable.hpp:134
void number(const int &number)
Set number of elements in column.
int m_inx5
Index of upper left node.
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
void append_table(const std::string &name, const std::string &unit)
Append table to response table.
Node array class.
Definition: GNodeArray.hpp:60
GCTAResponseTable(void)
Void constructor.
int m_ntables
Number of tables.
std::vector< std::string > m_units_hi
Upper boundaries units.
#define G_AXIS_RADIANS
GCTAResponseTable & operator=(const GCTAResponseTable &table)
Assignment operator.
void update(const double &arg) const
Update 1D cache.
#define G_UNIT
#define G_AXIS_LINEAR
double m_wgt6
Weight of lower left node.
#define G_ELEMENT_OPERATOR2
#define G_AXIS_LO_NAME
void axis_linear(const int &axis)
Set nodes for a linear axis.
double m_wgt2
Weight of lower left node.
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
std::vector< double > operator()(const double &arg) const
Linear interpolation operator for 1D tables.
#define G_INX_OPERATOR2
GFitsTableCol * append(const GFitsTableCol &column)
Append column to the table.
Definition: GFitsTable.hpp:147
void init_members(void)
Initialise response table members.
int m_inx2
Index of lower left node.
const int & tables(void) const
Return number of tables.
#define G_SCALE
#define G_AXIS_NODES
#define G_INX_OPERATOR3
int m_inx3
Index of upper right node.
void scale(const int &table, const double &scale)
Scale table.
void read_colnames(const GFitsTable &hdu)
Read column names from FITS HDU.
void read_axes(const GFitsTable &hdu)
Read axes definitions from FITS HDU.
Gammalib tools definition.
std::vector< std::vector< double > > m_axis_hi
Axes upper boundaries.
FITS table float column class interface definition.
#define G_AXIS_HI_NAME
void axis_radians(const int &axis)
Set nodes for a radians axis.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
#define G_AXIS_LOG10
#define G_APPEND_TABLE
virtual ~GCTAResponseTable(void)
Destructor.
void read(const GFitsTable &table)
Read response table from FITS table HDU.
#define G_OPERATOR1
#define G_OPERATOR2
#define G_APPEND_AXIS
int m_inx4
Index of lower right node.
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
std::vector< std::vector< double > > m_tables
Tables.
void write(GFitsTable &table) const
Write response table into FITS table HDU.
int axis(const std::string &name) const
Determine index of an axis.
#define G_TABLE
#define G_AXIS_LO
const std::string & unit(const int &table) const
Return table unit.
double m_wgt4
Weight of lower right node.
CTA response table class definition.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
bool has_table(const std::string &name) const
Check whether a table exists.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
#define G_READ_TABLES
void clear(void)
Clear response table.
Abstract interface for FITS table column.
std::vector< std::string > m_colname_table
Column names for table.
double m_wgt3
Weight of upper right node.
int m_inx7
Index of upper right node.
const int & elements(void) const
Return number of elements per table.
#define G_ELEMENT_OPERATOR1
GCTAResponseTable * clone(void) const
Clone response table.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
double m_wgt_left
Weight of left node.
GChatter
Definition: GTypemaps.hpp:33
void dim(const std::vector< int > &dim)
Set column dimension.
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
std::string print(const GChatter &chatter=NORMAL) const
Print response table information.
void free_members(void)
Delete response table members.
void axis_log10(const int &axis)
Set nodes for a logarithmic (base 10) axis.
#define G_READ_AXES
int table(const std::string &name) const
Determine index of table.
double m_wgt1
Weight of upper left node.
#define G_AXIS_LO_UNIT
void copy_members(const GCTAResponseTable &table)
Copy response table members.
double m_wgt8
Weight of lower right node.
#define G_OPERATOR3
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
std::vector< std::string > m_units_table
Parameter units.
double m_wgt7
Weight of upper right node.
#define G_INX_OPERATOR1
int axis_bins(const int &axis) const
Return number bins in an axis.
std::vector< std::vector< double > > m_axis_lo
Axes lower boundaries.
void append_axis(const std::vector< double > &axis_lo, const std::vector< double > &axis_hi, const std::string &name, const std::string &unit)
Append an axis to the response table.
const std::string & axis_lo_name(const int &axis) const
Return lower bin boundary FITS table column name for axis.
int m_inx_right
Index of right node.
std::vector< std::string > m_colname_hi
Column names for upper boundaries.
virtual double real(const int &row, const int &inx=0) const =0
int m_inx6
Index of lower left node.
std::vector< std::string > m_colname_lo
Column names for lower boundaries.
bool has_axis(const std::string &name) const
Check whether an axis exists.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
const std::string & axis_hi_name(const int &axis) const
Return upper bin boundary FITS table column name for axis.
#define G_AXIS_BINS
std::string string(const std::string &keyname) const
Return card value as string.
Definition: GFitsHDU.hpp:410
int m_naxes
Number of axes.
Exception handler interface definition.
std::vector< GNodeArray > m_axis_nodes
Axes node arrays.
void read_tables(const GFitsTable &hdu)
Read tables.
FITS binary table class definition.
FITS table float column.
CTA response table class.
#define G_READ_COLNAMES
int m_inx1
Index of upper left node.
const std::string & axis_lo_unit(const int &axis) const
Return lower bin boundary unit for axis.
#define G_AXIS_HI
const std::string & axis_hi_unit(const int &axis) const
Return upper bin boundary unit for axis.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
double m_wgt5
Weight of upper left node.
#define G_AXIS_HI_UNIT
int m_inx8
Index of lower right node.
GFitsHeaderCard & card(const int &cardno)
Return header card.
Definition: GFitsHDU.hpp:259
const int & axes(void) const
Return number of axes of the tables.
int m_nelements
Number of elements per table.
std::vector< std::string > m_units_lo
Lower boundaries units.
double m_wgt_right
Weight of right node.
Mathematical function definitions.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1205
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
std::string m_telescope
Telescope keyword.
#define G_AXIS