Review CTA effective area classes.

GCTAAeff2D The read() method now takes a reference instead of a pointer. GCTAAeffArf Rename read_arf() to read(). read() now takes a reference as argument instead of a pointer.

Review CTA effective area classes.
GCTAAeff2D The read() method now takes a reference instead of a pointer. GCTAAeffArf Rename read_arf() to read(). read() now takes a reference as argument instead of a pointer.
Jürgen Knödlseder
1 parent 45bcd2e8
Showing 14 changed files with 373 additions and 251 deletions
NEWS
inst/cta/include/GCTAAeff2D.hpp
inst/cta/include/GCTAAeffArf.hpp
inst/cta/include/GCTAAeffPerfTable.hpp
inst/cta/include/GCTAResponseTable.hpp
inst/cta/pyext/GCTAAeff2D.i
inst/cta/pyext/GCTAAeffArf.i
inst/cta/pyext/GCTAAeffPerfTable.i
inst/cta/pyext/GCTAResponseTable.i
inst/cta/src/GCTAAeff2D.cpp
inst/cta/src/GCTAAeffArf.cpp
inst/cta/src/GCTAAeffPerfTable.cpp
inst/cta/src/GCTAPsf2D.cpp
inst/cta/src/GCTAResponseTable.cpp
@@ -267,6 +267,13 @@ GWcsHPX =&gt; GHealpix
 GWcsRegistry
   Now holds objects of renamed type GWcs (see GHealpix).
+GCTAAeff2D
+  The read() method now takes a reference instead of a pointer.
+
+GCTAAeffArf
+  Rename read_arf() to read(). read() now takes a reference as
+  argument instead of a pointer.
+
 All GammaLib constants, defined in GTools.hpp, and functions are now
 put in the gammalib namespace to avoid name clashes. Mathematical
 constants and functions have been moved to a new model named GMath.
@@ -47,7 +47,7 @@ class GCTAAeff2D : public GCTAAeff {
 public:
     // Constructors and destructors
     GCTAAeff2D(void);
-    GCTAAeff2D(const std::string& filename);
+    explicit GCTAAeff2D(const std::string& filename);
     GCTAAeff2D(const GCTAAeff2D& cta);
     virtual ~GCTAAeff2D(void);
@@ -68,7 +68,7 @@ public:
     std::string print(const GChatter& chatter = NORMAL) const;
     // Methods
-    void read(const GFits* file);
+    void read(const GFits& file);
 private:
     // Methods
@@ -52,7 +52,7 @@ class GCTAAeffArf : public GCTAAeff {
 public:
     // Constructors and destructors
     GCTAAeffArf(void);
-    GCTAAeffArf(const std::string& filename);
+    explicit GCTAAeffArf(const std::string& filename);
     GCTAAeffArf(const GCTAAeffArf& cta);
     virtual ~GCTAAeffArf(void);
@@ -73,14 +73,14 @@ public:
     std::string  print(const GChatter& chatter = NORMAL) const;
     // Methods
-    int           size(void) const { return m_logE.size(); }
-    void          sigma(const double& sigma) { m_sigma=sigma; }
-    const double& sigma(void) const { return m_sigma; }
-    void          thetacut(const double& thetacut) { m_thetacut=thetacut; }
-    const double& thetacut(void) const { return m_thetacut; }
-    void          scale(const double& scale) { m_scale=scale; }
-    const double& scale(void) const { return m_scale; }
-    void          read_arf(const GFitsTable* hdu);
+    int           size(void) const;
+    void          sigma(const double& sigma);
+    const double& sigma(void) const;
+    void          thetacut(const double& thetacut);
+    const double& thetacut(void) const;
+    void          scale(const double& scale);
+    const double& scale(void) const;
+    void          read(const GFitsTable& hdu);
     void          remove_thetacut(const GCTAResponse& rsp);
 private:
@@ -98,4 +98,113 @@ private:
     double              m_scale;     //!< Scale for ARF
 };
+
+/***********************************************************************//**
+ * @brief Return filename
+ *
+ * @return Returns filename from which effective area was loaded.
+ ***************************************************************************/
+inline
+std::string GCTAAeffArf::filename(void) const
+{
+    return m_filename;
+}
+
+
+/***********************************************************************//**
+ * @brief Return number of node energies in response
+ *
+ * @return Number of node energies.
+ ***************************************************************************/
+inline
+int GCTAAeffArf::size(void) const
+{
+    return (m_logE.size());
+}
+
+
+/***********************************************************************//**
+ * @brief Set sigma for offset angle dependence
+ *
+ * @param[in] sigma Sigma for offset angle dependence.
+ *
+ * Sets the sigma parameter for the offset angle dependence of the effective
+ * area. If @p sigma is 0 (which is the default value), no offset angle
+ * dependency will be assumed.
+ ***************************************************************************/
+inline
+void GCTAAeffArf::sigma(const double& sigma)
+{
+    m_sigma = sigma;
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return sigma for offset angle dependence
+ *
+ * @return Sigma for offset angle dependence.
+ ***************************************************************************/
+inline
+const double& GCTAAeffArf::sigma(void) const
+{
+    return (m_sigma);
+}
+
+
+/***********************************************************************//**
+ * @brief Set theta cut angle
+ *
+ * @param[in] thetacut Set theta cut angle.
+ *
+ * Sets the theta cut angle which defines the energy independent cut that
+ * has been assumed in deriving the ARF values. If @p thetacut os 0 (which
+ * is the default avlue), not thetacut will be applied.
+ ***************************************************************************/
+inline
+void GCTAAeffArf::thetacut(const double& thetacut)
+{
+    m_thetacut = thetacut;
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return theta cut angle
+ *
+ * @return Theta cut angle.
+ ***************************************************************************/
+inline
+const double& GCTAAeffArf::thetacut(void) const
+{
+    return (m_thetacut);
+}
+
+
+/***********************************************************************//**
+ * @brief Set effective area scaling factor
+ *
+ * @param[in] scale Set effective area scaling factor.
+ *
+ * Sets the scaling factor for the effective area (by default is 1).
+ ***************************************************************************/
+inline
+void GCTAAeffArf::scale(const double& scale)
+{
+    m_scale = scale;
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return effective area scaling factor
+ *
+ * @return Effective area scaling factor.
+ ***************************************************************************/
+inline
+const double& GCTAAeffArf::scale(void) const
+{
+    return (m_scale);
+}
+
 #endif /* GCTAAEFFARF_HPP */
@@ -49,7 +49,7 @@ class GCTAAeffPerfTable : public GCTAAeff {
 public:
     // Constructors and destructors
     GCTAAeffPerfTable(void);
-    GCTAAeffPerfTable(const std::string& filename);
+    explicit GCTAAeffPerfTable(const std::string& filename);
     GCTAAeffPerfTable(const GCTAAeffPerfTable& cta);
     virtual ~GCTAAeffPerfTable(void);
@@ -70,9 +70,9 @@ public:
     std::string        print(const GChatter& chatter = NORMAL) const;
     // Methods
-    int           size(void) const { return m_logE.size(); }
-    void          sigma(const double& sigma) { m_sigma=sigma; }
-    const double& sigma(void) const { return m_sigma; }
+    int           size(void) const;
+    void          sigma(const double& sigma);
+    const double& sigma(void) const;
 private:
     // Methods
@@ -87,4 +87,57 @@ private:
     double              m_sigma;     //!< Sigma for offset angle computation (0=none)
 };
+
+/***********************************************************************//**
+ * @brief Return filename
+ *
+ * @return Returns filename from which effective area was loaded.
+ ***************************************************************************/
+inline
+std::string GCTAAeffPerfTable::filename(void) const
+{
+    return m_filename;
+}
+
+
+/***********************************************************************//**
+ * @brief Return number of node energies in response
+ *
+ * @return Number of node energies.
+ ***************************************************************************/
+inline
+int GCTAAeffPerfTable::size(void) const
+{
+    return (m_logE.size());
+}
+
+
+/***********************************************************************//**
+ * @brief Set sigma for offset angle dependence
+ *
+ * @param[in] sigma Sigma for offset angle dependence.
+ *
+ * Sets the sigma parameter for the offset angle dependence of the effective
+ * area. If @p sigma is 0, no offset angle dependency will be assumed. By
+ * default, @p sigma = 3.
+ ***************************************************************************/
+inline
+void GCTAAeffPerfTable::sigma(const double& sigma)
+{
+    m_sigma = sigma;
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Return sigma for offset angle dependence
+ *
+ * @return Sigma for offset angle dependence.
+ ***************************************************************************/
+inline
+const double& GCTAAeffPerfTable::sigma(void) const
+{
+    return (m_sigma);
+}
+
 #endif /* GCTAAEFFPERFTABLE_HPP */
@@ -49,11 +49,11 @@ public:
     // Constructors and destructors
     GCTAResponseTable(void);
     GCTAResponseTable(const GCTAResponseTable& table);
-    GCTAResponseTable(const GFitsTable* hdu);
+    explicit GCTAResponseTable(const GFitsTable& hdu);
     virtual ~GCTAResponseTable(void);
     // Operators
-    GCTAResponseTable&  operator= (const GCTAResponseTable & table);
+    GCTAResponseTable&  operator= (const GCTAResponseTable& table);
     std::vector<double> operator()(const double& arg) const;
     std::vector<double> operator()(const double& arg1, const double& arg2) const;
     double              operator()(const int& index, const double& arg) const;
@@ -73,7 +73,7 @@ public:
     void               axis_log10(const int& index);
     void               axis_radians(const int& index);
     void               scale(const int& index, const double& scale);
-    void               read(const GFitsTable* hdu);
+    void               read(const GFitsTable& hdu);
     void               write(GFitsTable* hdu) const;
     std::string        print(const GChatter& chatter = NORMAL) const;
@@ -82,9 +82,9 @@ private:
     void init_members(void);
     void copy_members(const GCTAResponseTable& table);
     void free_members(void);
-    void read_colnames(const GFitsTable* hdu);
-    void read_axes(const GFitsTable* hdu);
-    void read_pars(const GFitsTable* hdu);
+    void read_colnames(const GFitsTable& hdu);
+    void read_axes(const GFitsTable& hdu);
+    void read_pars(const GFitsTable& hdu);
     void update(const double& arg) const;
     void update(const double& arg1, const double& arg2) const;
@@ -40,7 +40,7 @@ class GCTAAeff2D : public GCTAAeff {
 public:
     // Constructors and destructors
     GCTAAeff2D(void);
-    GCTAAeff2D(const std::string& filename);
+    explicit GCTAAeff2D(const std::string& filename);
     GCTAAeff2D(const GCTAAeff2D& cta);
     virtual ~GCTAAeff2D(void);
@@ -59,7 +59,7 @@ public:
     std::string filename(void) const;
     // Methods
-    void read(const GFits* file);
+    void read(const GFits& file);
 };
 /***************************************************************************
  *                 GCTAAeffArf.i - CTA ARF effective area class            *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2012 by Juergen Knoedlseder                              *
+ *  copyright (C) 2012-2013 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -40,7 +40,7 @@ class GCTAAeffArf : public GCTAAeff {
 public:
     // Constructors and destructors
     GCTAAeffArf(void);
-    GCTAAeffArf(const std::string& filename);
+    explicit GCTAAeffArf(const std::string& filename);
     GCTAAeffArf(const GCTAAeffArf& cta);
     virtual ~GCTAAeffArf(void);
@@ -66,7 +66,8 @@ public:
     const double& thetacut(void) const;
     void          scale(const double& scale);
     const double& scale(void) const;
-    void          read_arf(const GFitsTable* hdu);
+    void          read(const GFitsTable& hdu);
+    void          remove_thetacut(const GCTAResponse& rsp);
 };
@@ -40,7 +40,7 @@ class GCTAAeffPerfTable : public GCTAAeff {
 public:
     // Constructors and destructors
     GCTAAeffPerfTable(void);
-    GCTAAeffPerfTable(const std::string& filename);
+    explicit GCTAAeffPerfTable(const std::string& filename);
     GCTAAeffPerfTable(const GCTAAeffPerfTable& cta);
     virtual ~GCTAAeffPerfTable(void);
@@ -46,7 +46,7 @@ public:
     // Constructors and destructors
     GCTAResponseTable(void);
     GCTAResponseTable(const GCTAResponseTable& table);
-    GCTAResponseTable(const GFitsTable* hdu);
+    explicit GCTAResponseTable(const GFitsTable& hdu);
     virtual ~GCTAResponseTable(void);
     // Interpolation operators
@@ -69,7 +69,7 @@ public:
     void               axis_log10(const int& index);
     void               axis_radians(const int& index);
     void               scale(const int& index, const double& scale);
-    void               read(const GFitsTable* hdu);
+    void               read(const GFitsTable& hdu);
     void               write(GFitsTable* hdu) const;
 };
@@ -243,7 +243,7 @@ void GCTAAeff2D::load(const std::string&amp; filename)
     GFits fits(filename);
     // Read effective area from file
-    read(&fits);
+    read(fits);
     // Close FITS file
     fits.close();
@@ -277,16 +277,16 @@ std::string GCTAAeff2D::filename(void) const
  * The data are stored in m_aeff which is of type GCTAResponseTable. The
  * energy axis will be set to log10, the offset angle axis to radians.
  ***************************************************************************/
-void GCTAAeff2D::read(const GFits* fits)
+void GCTAAeff2D::read(const GFits& fits)
 {
     // Clear response table
     m_aeff.clear();
-    // Get pointer to effective area HDU
-    const GFitsTable* hdu = fits->table("EFFECTIVE AREA");
+    // Get effective area table
+    const GFitsTable& table = *fits.table("EFFECTIVE AREA");
     // Read effective area table
-    m_aeff.read(hdu);
+    m_aeff.read(table);
     // Set energy axis to logarithmic scale
     m_aeff.axis_log10(0);
@@ -39,7 +39,6 @@
 #include "GCTAException.hpp"
 /* __ Method name definitions ____________________________________________ */
-#define G_LOAD                              "GCTAAeffArf::load(std::string&)"
 /* __ Macros _____________________________________________________________ */
@@ -134,7 +133,7 @@ GCTAAeffArf::~GCTAAeffArf(void)
  * @param[in] aeff Effective area.
  * @return Effective area.
  ***************************************************************************/
-GCTAAeffArf& GCTAAeffArf::operator= (const GCTAAeffArf& aeff)
+GCTAAeffArf& GCTAAeffArf::operator=(const GCTAAeffArf& aeff)
 {
     // Execute only if object is not identical
     if (this != &aeff) {
@@ -243,9 +242,6 @@ GCTAAeffArf* GCTAAeffArf::clone(void) const
  *
  * @param[in] filename Performance table file name.
  *
- * @exception GCTAExceptionHandler::file_open_error
- *            File could not be opened for read access.
- *
  * This method loads the effective area information from an ASCII
  * performance table.
  ***************************************************************************/
@@ -255,10 +251,10 @@ void GCTAAeffArf::load(const std::string&amp; filename)
     GFits file(filename);
     // Get ARF table
-    GFitsTable* table = file.table("SPECRESP");
+    const GFitsTable& table = *file.table("SPECRESP");
     // Read ARF
-    read_arf(table);
+    read(table);
     // Close ARF FITS file
     file.close();
@@ -272,21 +268,9 @@ void GCTAAeffArf::load(const std::string&amp; filename)
 /***********************************************************************//**
- * @brief Return filename
- *
- * @return Returns filename from which effective area was loaded
- ***************************************************************************/
-std::string GCTAAeffArf::filename(void) const
-{
-    // Return filename
-    return m_filename;
-}
-
-
-/***********************************************************************//**
  * @brief Read CTA ARF vector
  *
- * @param[in] hdu FITS table pointer.
+ * @param[in] hdu FITS table.
  *
  * This method reads a CTA ARF vector from the FITS HDU. Note that the
  * energies are converted to TeV and the effective area is converted to cm2.
@@ -298,16 +282,16 @@ std::string GCTAAeffArf::filename(void) const
  *       For appropriate theta angle assignment, we would need this
  *       information in the response header.
  ***************************************************************************/
-void GCTAAeffArf::read_arf(const GFitsTable* hdu)
+void GCTAAeffArf::read(const GFitsTable& hdu)
 {
     // Clear arrays
     m_logE.clear();
     m_aeff.clear();
     // Get pointers to table columns
-    const GFitsTableCol* energy_lo = (*hdu)["ENERG_LO"];
-    const GFitsTableCol* energy_hi = (*hdu)["ENERG_HI"];
-    const GFitsTableCol* specresp  = (*hdu)["SPECRESP"];
+    const GFitsTableCol* energy_lo = hdu["ENERG_LO"];
+    const GFitsTableCol* energy_hi = hdu["ENERG_HI"];
+    const GFitsTableCol* specresp  = hdu["SPECRESP"];
     // Determine unit conversion factors (default: TeV and cm^2)
     std::string u_energy_lo = gammalib::tolower(gammalib::strip_whitespace(energy_lo->unit()));
@@ -130,7 +130,7 @@ GCTAAeffPerfTable::~GCTAAeffPerfTable(void)
  * @param[in] aeff Effective area.
  * @return Effective area.
  ***************************************************************************/
-GCTAAeffPerfTable& GCTAAeffPerfTable::operator= (const GCTAAeffPerfTable& aeff)
+GCTAAeffPerfTable& GCTAAeffPerfTable::operator=(const GCTAAeffPerfTable& aeff)
 {
     // Execute only if object is not identical
     if (this != &aeff) {
@@ -303,18 +303,6 @@ void GCTAAeffPerfTable::load(const std::string&amp; filename)
 /***********************************************************************//**
- * @brief Return filename
- *
- * @return Returns filename from which effective area was loaded
- ***************************************************************************/
-std::string GCTAAeffPerfTable::filename(void) const
-{
-    // Return filename
-    return m_filename;
-}
-
-
-/***********************************************************************//**
  * @brief Print effective area information
  *
  * @param[in] chatter Chattiness (defaults to NORMAL).
@@ -259,7 +259,7 @@ void GCTAPsf2D::load(const std::string&amp; filename)
     GFits file(filename);
     // Get PSF table
-    GFitsTable* table = file.table("POINT SPREAD FUNCTION");
+    const GFitsTable& table = *file.table("POINT SPREAD FUNCTION");
     // Read PSF table
     m_psf.read(table);
@@ -119,7 +119,7 @@ GCTAResponseTable::GCTAResponseTable(const GCTAResponseTable&amp; table)
  * the object to a clean state (equivalent of GCTAResponseTable::clean
  * method).
  ***************************************************************************/
-GCTAResponseTable::GCTAResponseTable(const GFitsTable* hdu)
+GCTAResponseTable::GCTAResponseTable(const GFitsTable& hdu)
 {
     // Initialise class members
     init_members();
@@ -715,24 +715,19 @@ void GCTAResponseTable::scale(const int&amp; index, const double&amp; scale)
  * In case that the HDU table pointer is not valid (i.e. NULL), this method
  * clears the objects and does nothing else.
  ***************************************************************************/
-void GCTAResponseTable::read(const GFitsTable* hdu)
+void GCTAResponseTable::read(const GFitsTable& hdu)
 {
     // Clear instance
     clear();
-    // Continue only if HDU pointer is valid
-    if (hdu != NULL) {
+    // Read column names
+    read_colnames(hdu);
-        // Read column names
-        read_colnames(hdu);
+    // Read axes
+    read_axes(hdu);
-        // Read axes
-        read_axes(hdu);
-
-        // Read parameter cubes
-        read_pars(hdu);
-
-    } // endif: HDU pointer was valid
+    // Read parameter cubes
+    read_pars(hdu);
     // Return
     return;
@@ -923,7 +918,7 @@ void GCTAResponseTable::free_members(void)
 /***********************************************************************//**
  * @brief Read column names from FITS HDU
  *
- * @param[in] hdu Response table HDU pointer.
+ * @param[in] hdu Response table HDU.
  *
  * @exception GCTAException::bad_rsp_table_format
  *            Bad response table format encountered in FITS HDU.
@@ -949,7 +944,7 @@ void GCTAResponseTable::free_members(void)
  *
  * @todo Implement exceptions for invalid HDU format
  ***************************************************************************/
-void GCTAResponseTable::read_colnames(const GFitsTable* hdu)
+void GCTAResponseTable::read_colnames(const GFitsTable& hdu)
 {
     // Clear column name arrays
     m_naxes = 0;
@@ -958,103 +953,98 @@ void GCTAResponseTable::read_colnames(const GFitsTable* hdu)
     m_colname_hi.clear();
     m_colname_par.clear();
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Initialise search mode. There are three search modes:
-        // 0 - we're looking for the next axis by searching for a column
-        //     terminating with "_LO"
-        // 1 - we're looking for the upper boundary of an axis, terminating
-        //     with "_HI"
-        // 2 - we're looking for a parameter column
-        int         mode = 0;
-        std::string lo_column;
-        std::string next_column;
-
-        // Extract column names for all axes
-        for (int i = 0; i < hdu->ncols(); ++i) {
-
-            // Get column name
-            std::string colname = (*hdu)[i]->name();
-
-            // If we search for a "_LO" column, check if we have one. If one
-            // is found, change the search mode to 1 and set the expected name
-            // for the "_HI" column. If none is found, change the search
-            // mode to 2 since from now on we should only have parameter
-            // columns.
-            if (mode == 0) {
-                size_t pos = colname.rfind("_LO");
-                if (pos != std::string::npos) {
-                    mode        = 1;
-                    lo_column   = colname;
-                    next_column = colname.substr(0, pos) + "_HI";
-                }
-                else {
-                    if (colname.rfind("_HI") != std::string::npos) {
-                        std::string message = "Column '" +
-                                              colname +
-                                              "' encountered without"
-                                              " preceeding '_LO' column.";
-                        throw GCTAException::bad_rsp_table_format(G_READ_COLNAMES,
-                                                                  message);
-                    }
-                    else {
-                        mode = 2;
-                        m_colname_par.push_back(colname);
-                    }
-                }
-            }
-
-            // If we search for a "_HI" column, check if we have the
-            // expected column name. If this is the case, switch back to
-            // search mode 0 to get the next "_LO" column. Otherwise
-            // throw an exception.
-            else if (mode == 1) {
-                if (colname == next_column) {
-                    mode = 0;
-                    m_colname_lo.push_back(lo_column);
-                    m_colname_hi.push_back(next_column);
-                }
-                else {
-                    std::string message = "Expected column '" +
-                                          next_column +
-                                          "' not found. '_HI' columns have"
-                                          " to be placed immediately after"
-                                          " corresponding '_LO' columns.";
-                    throw GCTAException::bad_rsp_table_format(G_READ_COLNAMES,
-                                                              message);
-                }
+    // Initialise search mode. There are three search modes:
+    // 0 - we're looking for the next axis by searching for a column
+    //     terminating with "_LO"
+    // 1 - we're looking for the upper boundary of an axis, terminating
+    //     with "_HI"
+    // 2 - we're looking for a parameter column
+    int         mode = 0;
+    std::string lo_column;
+    std::string next_column;
+
+    // Extract column names for all axes
+    for (int i = 0; i < hdu.ncols(); ++i) {
+
+        // Get column name
+        std::string colname = hdu[i]->name();
+
+        // If we search for a "_LO" column, check if we have one. If one
+        // is found, change the search mode to 1 and set the expected name
+        // for the "_HI" column. If none is found, change the search
+        // mode to 2 since from now on we should only have parameter
+        // columns.
+        if (mode == 0) {
+            size_t pos = colname.rfind("_LO");
+            if (pos != std::string::npos) {
+                mode        = 1;
+                lo_column   = colname;
+                next_column = colname.substr(0, pos) + "_HI";
             }
-
-            // If we search for a parameter column, make sure that we have
-            // neither a "_LO" nor a "_HI" column.
             else {
-                if (colname.rfind("_LO") != std::string::npos) {
+                if (colname.rfind("_HI") != std::string::npos) {
                     std::string message = "Column '" +
                                           colname +
-                                          "' found. All '_LO' columns have to"
-                                          " be placed before the parameter"
-                                          " columns.";
-                    throw GCTAException::bad_rsp_table_format(G_READ_COLNAMES,
-                                                              message);
-                }
-                else if (colname.rfind("_HI") != std::string::npos) {
-                    std::string message = "Column '" +
-                                          colname +
-                                          "' found. All '_HI' columns have to"
-                                          " be placed before the parameter"
-                                          " columns.";
+                                          "' encountered without"
+                                          " preceeding '_LO' column.";
                     throw GCTAException::bad_rsp_table_format(G_READ_COLNAMES,
                                                               message);
                 }
                 else {
+                    mode = 2;
                     m_colname_par.push_back(colname);
                 }
             }
+        }
-        } // endfor: looped over all column names
+        // If we search for a "_HI" column, check if we have the
+        // expected column name. If this is the case, switch back to
+        // search mode 0 to get the next "_LO" column. Otherwise
+        // throw an exception.
+        else if (mode == 1) {
+            if (colname == next_column) {
+                mode = 0;
+                m_colname_lo.push_back(lo_column);
+                m_colname_hi.push_back(next_column);
+            }
+            else {
+                std::string message = "Expected column '" +
+                                      next_column +
+                                      "' not found. '_HI' columns have"
+                                      " to be placed immediately after"
+                                      " corresponding '_LO' columns.";
+                throw GCTAException::bad_rsp_table_format(G_READ_COLNAMES,
+                                                          message);
+            }
+        }
-    } // endif: HDU was valid
+        // If we search for a parameter column, make sure that we have
+        // neither a "_LO" nor a "_HI" column.
+        else {
+            if (colname.rfind("_LO") != std::string::npos) {
+                std::string message = "Column '" +
+                                      colname +
+                                      "' found. All '_LO' columns have to"
+                                      " be placed before the parameter"
+                                      " columns.";
+                throw GCTAException::bad_rsp_table_format(G_READ_COLNAMES,
+                                                          message);
+            }
+            else if (colname.rfind("_HI") != std::string::npos) {
+                std::string message = "Column '" +
+                                      colname +
+                                      "' found. All '_HI' columns have to"
+                                      " be placed before the parameter"
+                                      " columns.";
+                throw GCTAException::bad_rsp_table_format(G_READ_COLNAMES,
+                                                          message);
+            }
+            else {
+                m_colname_par.push_back(colname);
+            }
+        }
+
+    } // endfor: looped over all column names
     // Store number of axes
     m_naxes = m_colname_lo.size();
@@ -1070,7 +1060,7 @@ void GCTAResponseTable::read_colnames(const GFitsTable* hdu)
 /***********************************************************************//**
  * @brief Read axes definitions from FITS HDU
  *
- * @param[in] hdu Response table HDU pointer.
+ * @param[in] hdu Response table HDU.
  *
  * @exception GCTAException::bad_rsp_table_format
  *            Incompatible axes columns found.
@@ -1094,58 +1084,53 @@ void GCTAResponseTable::read_colnames(const GFitsTable* hdu)
  * In case that the HDU pointer is not valid (NULL), this method clears the
  * axes boundaries and does nothing else.
  ***************************************************************************/
-void GCTAResponseTable::read_axes(const GFitsTable* hdu)
+void GCTAResponseTable::read_axes(const GFitsTable& hdu)
 {
     // Clear axes arrays
     m_axis_lo.clear();
     m_axis_hi.clear();
     m_axis_nodes.clear();
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Loop over all dimensions
-        for (int i = 0; i < axes(); ++i) {
-
-            // Get pointers to table columns
-            const GFitsTableCol* col_lo = (*hdu)[m_colname_lo[i]];
-            const GFitsTableCol* col_hi = (*hdu)[m_colname_hi[i]];
-
-            // Extract number of bins. Make sure that both columns have the
-            // same number of bins
-            int num = col_lo->number();
-            if (num != col_hi->number()) {
-                std::string message = "Incompatible number of bins found"
-                                      " for columns '"+m_colname_lo[i]+"'"+
-                                      " ("+gammalib::str(num)+") and '"+
-                                      m_colname_hi[i]+"' ("+
-                                      gammalib::str(col_hi->number())+".";
-                throw GCTAException::bad_rsp_table_format(G_READ_AXES,
-                                                          message);
-            }
-
-            // Initialise axis and node arrays
-            std::vector<double> axis_lo(num);
-            std::vector<double> axis_hi(num);
-            std::vector<double> axis_nodes(num);
+    // Loop over all dimensions
+    for (int i = 0; i < axes(); ++i) {
+
+        // Get pointers to table columns
+        const GFitsTableCol* col_lo = hdu[m_colname_lo[i]];
+        const GFitsTableCol* col_hi = hdu[m_colname_hi[i]];
+
+        // Extract number of bins. Make sure that both columns have the
+        // same number of bins
+        int num = col_lo->number();
+        if (num != col_hi->number()) {
+            std::string message = "Incompatible number of bins found"
+                                  " for columns '"+m_colname_lo[i]+"'"+
+                                  " ("+gammalib::str(num)+") and '"+
+                                  m_colname_hi[i]+"' ("+
+                                  gammalib::str(col_hi->number())+".";
+            throw GCTAException::bad_rsp_table_format(G_READ_AXES,
+                                                      message);
+        }
-            // Copy axis information into arrays
-            for (int k = 0; k < num; ++k) {
-                axis_lo[k]    = col_lo->real(0,k);
-                axis_hi[k]    = col_hi->real(0,k);
-                axis_nodes[k] = 0.5*(axis_lo[k] + axis_hi[k]);
-            }
+        // Initialise axis and node arrays
+        std::vector<double> axis_lo(num);
+        std::vector<double> axis_hi(num);
+        std::vector<double> axis_nodes(num);
-            // Push axis array on storage
-            m_axis_lo.push_back(axis_lo);
-            m_axis_hi.push_back(axis_hi);
+        // Copy axis information into arrays
+        for (int k = 0; k < num; ++k) {
+            axis_lo[k]    = col_lo->real(0,k);
+            axis_hi[k]    = col_hi->real(0,k);
+            axis_nodes[k] = 0.5*(axis_lo[k] + axis_hi[k]);
+        }
-            // Create node array
-            m_axis_nodes.push_back(GNodeArray(axis_nodes));
+        // Push axis array on storage
+        m_axis_lo.push_back(axis_lo);
+        m_axis_hi.push_back(axis_hi);
-        } // endfor: looped over all dimensions
+        // Create node array
+        m_axis_nodes.push_back(GNodeArray(axis_nodes));
-    } // endif: HDU was valid
+    } // endfor: looped over all dimensions
     // Return
     return;
@@ -1155,7 +1140,7 @@ void GCTAResponseTable::read_axes(const GFitsTable* hdu)
 /***********************************************************************//**
  * @brief Read parameter cubes
  *
- * @param[in] hdu Response table HDU pointer.
+ * @param[in] hdu Response table HDU.
  *
  * @exception GCTAException::bad_rsp_table_format
  *            Parameter vector of bad size encountered.
@@ -1171,52 +1156,47 @@ void GCTAResponseTable::read_axes(const GFitsTable* hdu)
  * In case that the HDU pointer is not valid (NULL), this method clears the
  * axes boundaries and does nothing else.
  ***************************************************************************/
-void GCTAResponseTable::read_pars(const GFitsTable* hdu)
+void GCTAResponseTable::read_pars(const GFitsTable& hdu)
 {
     // Clear parameter cubes
     m_pars.clear();
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Compute expected cube size
-        m_nelements = axis(0);
-        for (int i = 1; i < axes(); ++i) {
-            m_nelements *= axis(i);
-        }
+    // Compute expected cube size
+    m_nelements = axis(0);
+    for (int i = 1; i < axes(); ++i) {
+        m_nelements *= axis(i);
+    }
-        // Loop over all parameter cubes
-        for (int i = 0; i < size(); ++i) {
-
-            // Get pointer to table column
-            const GFitsTableCol* col = (*hdu)[m_colname_par[i]];
-
-            // Extract number of bins. Verify that the number of bins
-            // corresponds to the expectation.
-            int num = col->number();
-            if (num != m_nelements) {
-                std::string message = "Parameter vector '"+m_colname_par[i]+
-                                      "' has wrong size "+gammalib::str(num)+
-                                      " (expected"+
-                                      gammalib::str(m_nelements)+").";
-                throw GCTAException::bad_rsp_table_format(G_READ_PARS,
-                                                          message);
-            }
-
-            // Initialise parameter cube
-            std::vector<double> pars(num);
+    // Loop over all parameter cubes
+    for (int i = 0; i < size(); ++i) {
+
+        // Get pointer to table column
+        const GFitsTableCol* col = hdu[m_colname_par[i]];
+
+        // Extract number of bins. Verify that the number of bins
+        // corresponds to the expectation.
+        int num = col->number();
+        if (num != m_nelements) {
+            std::string message = "Parameter vector '"+m_colname_par[i]+
+                                  "' has wrong size "+gammalib::str(num)+
+                                  " (expected"+
+                                  gammalib::str(m_nelements)+").";
+            throw GCTAException::bad_rsp_table_format(G_READ_PARS,
+                                                      message);
+        }
-            // Copy parameter values
-            for (int k = 0; k < num; ++k) {
-                pars[k] = col->real(0,k);
-            }
+        // Initialise parameter cube
+        std::vector<double> pars(num);
-            // Push cube into storage
-            m_pars.push_back(pars);
+        // Copy parameter values
+        for (int k = 0; k < num; ++k) {
+            pars[k] = col->real(0,k);
+        }
-        } // endfor: looped over all parameter cubes
+        // Push cube into storage
+        m_pars.push_back(pars);
-    } // endif: HDU was valid
+    } // endfor: looped over all parameter cubes
     // Return
     return;