The FITS read() and write() methods now take reference instead of pointers.

Jürgen Knödlseder
1 parent 3e38d75b
Showing 32 changed files with 644 additions and 641 deletions
NEWS
include/GArf.hpp
include/GEnergies.hpp
include/GHealpix.hpp
include/GPha.hpp
include/GPhotons.hpp
include/GRmf.hpp
include/GSkyProjection.hpp
include/GSkymap.hpp
include/GTimes.hpp
include/GWcs.hpp
inst/com/src/GCOMEventCube.cpp
inst/com/src/GCOMObservation.cpp
inst/cta/pyext/GCTAInstDir.i
inst/cta/src/GCTAEventCube.cpp
inst/lat/src/GLATEventCube.cpp
inst/lat/src/GLATLtCubeMap.cpp
pyext/GArf.i
pyext/GEnergies.i
pyext/GHealpix.i
@@ -180,6 +180,13 @@ GPars
 GPhoton
   Add mcid argument to Photon constructor (defaults to -1).
+GSkymap
+  Rename xy2dir() to pix2dir(), dir2xy() to dir2pix(), xy2pix() to
+  pix2inx() and pix2xy() to inx2pix(). The nmaps argument in the
+  sky map constructor is now passed by reference. The isinmap()
+  method has been renamed to contains(). The read() and write()
+  methods now take a reference instead of a pointer.
+
 GSkyPixel
   The class now supports 1D and 2D pixels. Type conversion methods
   have been added that allow assigning integer and double precision
@@ -222,16 +229,19 @@ GWcs =&gt; GSkyProjection
   general than the World Coordinate Systems. A size() method has
   been added that returns the dimension of the projection (1D or
   2D). The integer index methods have been removed as GSkyPixel
-  now supports 1D and 2D pixels.
+  now supports 1D and 2D pixels. The read() and write() methods
+  now take a reference instead of a pointer.
 GWcslib => GWcs
   The GWcslib class has been renamed to GWcs. The class now conforms
-  to the GSkyProjection base class interface.
+  to the GSkyProjection base class interface. The read() and write()
+  methods now take a reference instead of a pointer.
 GWcsHPX => GHealpix
   The GWcsHPX class has been renamed to GHealpix as the Healpix
   projection is not a World Coordinate System. GHealpix is no
-  longer registered in GWcsRegistry.
+  longer registered in GWcsRegistry. The read() and write()
+  methods now take a reference instead of a pointer.
 GWcsRegistry
   Now holds objects of renamed type GWcs (see GHealpix).
@@ -63,8 +63,9 @@ public:
     const double&      at(const int& index) const;
     const GEbounds&    ebounds(void) const;
     void               load(const std::string& filename);
-    void               save(const std::string& filename, const bool& clobber = false) const;
-    void               read(const GFitsTable* hdu);
+    void               save(const std::string& filename,
+                            const bool& clobber = false) const;
+    void               read(const GFitsTable& table);
     void               write(GFits& fits) const;
     const std::string& filename(void) const;
     std::string        print(const GChatter& chatter = NORMAL) const;
@@ -75,8 +75,8 @@ public:
                         const std::string& extname = "ENERGIES");
     void           save(const std::string& filename, bool clobber,
                         const std::string& extname = "ENERGIES") const;
-    void           read(const GFitsTable* hdu);
-    void           write(GFits* file,
+    void           read(const GFitsTable& table);
+    void           write(GFits& file,
                          const std::string& extname = "ENERGIES") const;
     std::string    print(const GChatter& chatter = NORMAL) const;
@@ -51,7 +51,7 @@ public:
     explicit GHealpix(const int&         nside,
                       const std::string& ordering = "NESTED",
                       const std::string& coordsys = "GAL");
-    explicit GHealpix(const GFitsHDU* hdu);
+    explicit GHealpix(const GFitsHDU& hdu);
     GHealpix(const GHealpix& wcs);
     virtual ~GHealpix(void);
@@ -64,8 +64,8 @@ public:
     virtual int         size(void) const;
     virtual std::string code(void) const;
     virtual std::string name(void) const;
-    virtual void        read(const GFitsHDU* hdu);
-    virtual void        write(GFitsHDU* hdu) const;
+    virtual void        read(const GFitsHDU& hdu);
+    virtual void        write(GFitsHDU& hdu) const;
     virtual double      omega(const GSkyPixel& pixel) const;
     virtual GSkyDir     pix2dir(const GSkyPixel& pixel) const;
     virtual GSkyPixel   dir2pix(const GSkyDir& dir) const;
@@ -80,8 +80,9 @@ public:
     const double&      outflow(void) const;
     void               fill(const GEnergy& energy, const double& value = 1.0);
     void               load(const std::string& filename);
-    void               save(const std::string& filename, const bool& clobber = false) const;
-    void               read(const GFitsTable* hdu);
+    void               save(const std::string& filename,
+                            const bool& clobber = false) const;
+    void               read(const GFitsTable& table);
     void               write(GFits& fits) const;
     const std::string& filename(void) const;
     std::string        print(const GChatter& chatter = NORMAL) const;
@@ -58,8 +58,8 @@ public:
     // Methods
     void           clear(void);
     GPhotons*      clone(void) const;
-    int            size(void) const { return m_photons.size(); }
-    bool           isempty(void) const { return m_photons.empty(); }
+    int            size(void) const;
+    bool           isempty(void) const;
     void           append(const GPhoton& photon);
     void           insert(const int& index, const GPhoton& photon);
     void           remove(const int& index);
@@ -77,4 +77,28 @@ protected:
     std::vector<GPhoton> m_photons;  //!< List of photons
 };
+
+/***********************************************************************//**
+ * @brief Return number of photons
+ *
+ * @return Number of photons.
+ ***************************************************************************/
+inline
+int GPhotons::size(void) const
+{
+    return m_photons.size();
+}
+
+
+/***********************************************************************//**
+ * @brief Signal if there are no photons
+ *
+ * @return True if there are no photons.
+ ***************************************************************************/
+inline
+bool GPhotons::isempty(void) const
+{
+    return m_photons.empty();
+}
+
 #endif /* GPHOTONS_HPP */
@@ -69,8 +69,9 @@ public:
     const GEbounds&    etrue(void) const;
     const GEbounds&    emeasured(void) const;
     void               load(const std::string& filename);
-    void               save(const std::string& filename, const bool& clobber = false) const;
-    void               read(const GFitsTable* hdu);
+    void               save(const std::string& filename,
+                            const bool& clobber = false) const;
+    void               read(const GFitsTable& table);
     void               write(GFits& fits) const;
     const std::string& filename(void) const;
     std::string        print(const GChatter& chatter = NORMAL) const;
@@ -65,8 +65,8 @@ public:
     virtual int             size(void) const = 0;
     virtual std::string     code(void) const = 0;
     virtual std::string     name(void) const = 0;
-    virtual void            read(const GFitsHDU* hdu) = 0;
-    virtual void            write(GFitsHDU* hdu) const = 0;
+    virtual void            read(const GFitsHDU& hdu) = 0;
+    virtual void            write(GFitsHDU& hdu) const = 0;
     virtual double          omega(const GSkyPixel& pixel) const = 0;
     virtual GSkyDir         pix2dir(const GSkyPixel& pixel) const = 0;
     virtual GSkyPixel       dir2pix(const GSkyDir& dir) const = 0;
@@ -128,8 +128,8 @@ public:
     const double*         pixels(void) const;
     void                  load(const std::string& filename);
     void                  save(const std::string& filename, bool clobber = false) const;
-    void                  read(const GFitsHDU* hdu);
-    void                  write(GFits* file) const;
+    void                  read(const GFitsHDU& hdu);
+    void                  write(GFits& file) const;
     std::string           print(const GChatter& chatter = NORMAL) const;
 private:
@@ -143,9 +143,9 @@ private:
                               const double& crpix1, const double& crpix2,
                               const double& cdelt1, const double& cdelt2,
                               const GMatrix& cd, const GVector& pv2);
-    void              read_healpix(const GFitsTable* hdu);
-    void              read_wcs(const GFitsImage* hdu);
-    void              alloc_wcs(const GFitsImage* hdu);
+    void              read_healpix(const GFitsTable& table);
+    void              read_wcs(const GFitsImage& image);
+    void              alloc_wcs(const GFitsImage& image);
     GFitsBinTable*    create_healpix_hdu(void) const;
     GFitsImageDouble* create_wcs_hdu(void) const;
@@ -58,8 +58,8 @@ public:
     // Methods
     void        clear(void);
     GTimes*     clone(void) const;
-    int         size(void) const { return m_times.size(); }
-    bool        isempty(void) const { return m_times.empty(); }
+    int         size(void) const;
+    bool        isempty(void) const;
     void        append(const GTime& time);
     void        insert(const int& index, const GTime& time);
     void        remove(const int& index);
@@ -77,4 +77,28 @@ protected:
     std::vector<GTime> m_times;  //!< List of times
 };
+
+/***********************************************************************//**
+ * @brief Return number of times
+ *
+ * @return Number of times.
+ ***************************************************************************/
+inline
+int GTimes::size(void) const
+{
+    return m_times.size();
+}
+
+
+/***********************************************************************//**
+ * @brief Signal if there are no times
+ *
+ * @return True if there are no times.
+ ***************************************************************************/
+inline
+bool GTimes::isempty(void) const
+{
+    return m_times.empty();
+}
+
 #endif /* GTIMES_HPP */
@@ -52,7 +52,7 @@ public:
                   const double& crval1, const double& crval2,
                   const double& crpix1, const double& crpix2,
                   const double& cdelt1, const double& cdelt2);
-    explicit GWcs(const GFitsHDU* hdu);
+    explicit GWcs(const GFitsHDU& hdu);
     GWcs(const GWcs& wcs);
     virtual ~GWcs(void);
@@ -68,8 +68,8 @@ public:
     // Implemented virtual base class methods
     virtual int         size(void) const;
-    virtual void        read(const GFitsHDU* hdu);
-    virtual void        write(GFitsHDU* hdu) const;
+    virtual void        read(const GFitsHDU& hdu);
+    virtual void        write(GFitsHDU& hdu) const;
     virtual double      omega(const GSkyPixel& pixel) const;
     virtual GSkyDir     pix2dir(const GSkyPixel& pixel) const;
     virtual GSkyPixel   dir2pix(const GSkyDir& dir) const;
@@ -417,11 +417,11 @@ void GCOMEventCube::read(const GFits&amp; fits)
     // Clear object
     clear();
-    // Get HDU
-    const GFitsImage* hdu = fits.image("Primary");
+    // Get image
+    const GFitsImage& image = *fits.image("Primary");
     // Load counts map as sky map
-    m_map.read(hdu);
+    m_map.read(image);
     // Correct WCS projection (HEASARC data format kluge)
     com_wcs_mer2car(m_map);
@@ -430,9 +430,9 @@ void GCOMEventCube::read(const GFits&amp; fits)
     set_scatter_directions();
     // Extract Phibar minimum and step size
-    double crval3 = hdu->real("CRVAL3");
-    double crpix3 = hdu->real("CRPIX3");
-    double dphi   = hdu->real("CDELT3");
+    double crval3 = image.real("CRVAL3");
+    double crpix3 = image.real("CRPIX3");
+    double dphi   = image.real("CDELT3");
     double phimin = crval3 + (crpix3-1.0) * dphi;
     // Set scatter angles
@@ -441,8 +441,8 @@ void GCOMEventCube::read(const GFits&amp; fits)
     // Extract energy range
     GEnergy emin;
     GEnergy emax;
-    emin.MeV(hdu->real("E_MIN"));
-    emax.MeV(hdu->real("E_MAX"));
+    emin.MeV(image.real("E_MIN"));
+    emax.MeV(image.real("E_MAX"));
     // Set energy boundaries
     m_ebounds.clear();
@@ -454,8 +454,8 @@ void GCOMEventCube::read(const GFits&amp; fits)
     // Extract time range
     GTime  tstart;
     GTime  tstop;
-    tstart.jd(hdu->real("TSTART"));
-    tstop.jd(hdu->real("TSTOP"));
+    tstart.jd(image.real("TSTART"));
+    tstop.jd(image.real("TSTOP"));
     // Set GTIs
     m_gti.clear();
@@ -479,7 +479,7 @@ void GCOMEventCube::read(const GFits&amp; fits)
 void GCOMEventCube::write(GFits& file) const
 {
     // Write cube
-    m_map.write(&file);
+    m_map.write(file);
     // Write energy boundaries
@@ -734,11 +734,11 @@ void GCOMObservation::load_drb(const std::string&amp; drbname)
     // Open FITS file
     GFits file(drbname);
-    // Get HDU
-    GFitsImage* hdu = file.image("Primary");
+    // Get image
+    const GFitsImage& image = *file.image("Primary");
     // Load background model as sky map
-    m_drb.read(hdu);
+    m_drb.read(image);
     // Correct WCS projection (HEASARC data format kluge)
     com_wcs_mer2car(m_drb);
@@ -776,11 +776,11 @@ void GCOMObservation::load_drg(const std::string&amp; drgname)
     // Open FITS file
     GFits file(drgname);
-    // Get HDU
-    GFitsImage* hdu = file.image("Primary");
+    // Get image
+    const GFitsImage& image = *file.image("Primary");
     // Load geometry factors as sky map
-    m_drg.read(hdu);
+    m_drg.read(image);
     // Correct WCS projection (HEASARC data format kluge)
     com_wcs_mer2car(m_drg);
@@ -815,10 +815,10 @@ void GCOMObservation::load_drx(const std::string&amp; drxname)
     GFits file(drxname);
     // Get HDU
-    GFitsImage* hdu = file.image("Primary");
+    const GFitsImage& image = *file.image("Primary");
     // Load exposure map as sky map
-    m_drx.read(hdu);
+    m_drx.read(image);
     // Correct WCS projection (HEASARC data format kluge)
     com_wcs_mer2car(m_drx);
@@ -48,9 +48,8 @@ public:
     GCTAInstDir* clone(void) const;
     // Other methods
-    void           dir(const GSkyDir& dir);
-    GSkyDir&       dir(void);
-    const GSkyDir& dir(void) const;
+    void     dir(const GSkyDir& dir);
+    GSkyDir& dir(void);
 };
@@ -396,7 +396,7 @@ void GCTAEventCube::read(const GFits&amp; fits)
 void GCTAEventCube::write(GFits& fits) const
 {
     // Write cube
-    m_map.write(&fits);
+    m_map.write(fits);
     // Write energy boundaries
     ebounds().write(fits);
@@ -574,7 +574,7 @@ void GCTAEventCube::read_cntmap(const GFitsImage* hdu)
     if (hdu != NULL) {
         // Load counts map as sky map
-        m_map.read(hdu);
+        m_map.read(*hdu);
         // Set sky directions
         set_directions();
@@ -695,7 +695,7 @@ void GLATEventCube::read_cntmap(const GFitsImage* hdu)
     if (hdu != NULL) {
         // Load counts map as sky map
-        m_map.read(hdu);
+        m_map.read(*hdu);
         // Set sky directions
         set_directions();
@@ -727,7 +727,7 @@ void GLATEventCube::read_srcmap(const GFitsImage* hdu)
         GSkymap* map = new GSkymap;
         // Read skymap
-        map->read(hdu);
+        map->read(*hdu);
         // Check that source map sky projection is consistent with counts
         // map sky projection
@@ -372,7 +372,7 @@ void GLATLtCubeMap::read(const GFitsTable* hdu)
     clear();
     // Load skymap
-    m_map.read(hdu);
+    m_map.read(*hdu);
     // Set costheta binning scheme
     std::string scheme = gammalib::strip_whitespace(gammalib::toupper(hdu->string("THETABIN")));
@@ -51,8 +51,9 @@ public:
     double&         at(const int& index);
     const GEbounds& ebounds(void) const;
     void            load(const std::string& filename);
-    void            save(const std::string& filename, const bool& clobber = false) const;
-    void            read(const GFitsTable* hdu);
+    void            save(const std::string& filename,
+                         const bool& clobber = false) const;
+    void            read(const GFitsTable& table);
     void            write(GFits& fits) const;
 };
@@ -61,8 +61,8 @@ public:
                     const std::string& extname = "ENERGIES");
     void       save(const std::string& filename, bool clobber,
                     const std::string& extname = "ENERGIES") const;
-    void       read(const GFitsTable* hdu);
-    void       write(GFits* file, const std::string& extname = "ENERGIES") const;
+    void       read(const GFitsTable& table);
+    void       write(GFits& file, const std::string& extname = "ENERGIES") const;
 };
@@ -40,7 +40,7 @@ public:
     GHealpix(void);
     explicit GHealpix(const int& nside, const std::string& ordering = "NESTED",
                       const std::string& coordsys = "GAL");
-    explicit GHealpix(const GFitsHDU* hdu);
+    explicit GHealpix(const GFitsHDU& hdu);
     GHealpix(const GHealpix& wcs);
     virtual ~GHealpix(void);
@@ -50,8 +50,8 @@ public:
     virtual int         size(void) const;
     virtual std::string code(void) const;
     virtual std::string name(void) const;
-    virtual void        read(const GFitsHDU* hdu);
-    virtual void        write(GFitsHDU* hdu) const;
+    virtual void        read(const GFitsHDU& hdu);
+    virtual void        write(GFitsHDU& hdu) const;
     virtual double      omega(const GSkyPixel& pixel) const;
     virtual GSkyDir     pix2dir(const GSkyPixel& pixel) const;
     virtual GSkyPixel   dir2pix(const GSkyDir& dir) const;
@@ -57,8 +57,9 @@ public:
     const double&      outflow(void) const;
     void               fill(const GEnergy& energy, const double& value = 1.0);
     void               load(const std::string& filename);
-    void               save(const std::string& filename, const bool& clobber = false) const;
-    void               read(const GFitsTable* hdu);
+    void               save(const std::string& filename,
+                            const bool& clobber = false) const;
+    void               read(const GFitsTable& table);
     void               write(GFits& fits) const;
     const std::string& filename(void) const;
 };
@@ -57,8 +57,9 @@ public:
     const GEbounds& etrue(void) const;
     const GEbounds& emeasured(void) const;
     void            load(const std::string& filename);
-    void            save(const std::string& filename, const bool& clobber = false) const;
-    void            read(const GFitsTable* hdu);
+    void            save(const std::string& filename,
+                         const bool& clobber = false) const;
+    void            read(const GFitsTable& table);
     void            write(GFits& fits) const;
 };
@@ -51,8 +51,8 @@ public:
     virtual int             size(void) const = 0;
     virtual std::string     code(void) const = 0;
     virtual std::string     name(void) const = 0;
-    virtual void            read(const GFitsHDU* hdu) = 0;
-    virtual void            write(GFitsHDU* hdu) const = 0;
+    virtual void            read(const GFitsHDU& hdu) = 0;
+    virtual void            write(GFitsHDU& hdu) const = 0;
     virtual double          omega(const GSkyPixel& pixel) const = 0;
     virtual GSkyDir         pix2dir(const GSkyPixel& pixel) const = 0;
     virtual GSkyPixel       dir2pix(const GSkyDir& dir) const = 0;
@@ -149,8 +149,8 @@ public:
     const double*         pixels(void) const;
     void                  load(const std::string& filename);
     void                  save(const std::string& filename, bool clobber = false) const;
-    void                  read(const GFitsHDU* hdu);
-    void                  write(GFits* file) const;
+    void                  read(const GFitsHDU& hdu);
+    void                  write(GFits& file) const;
 };
@@ -42,7 +42,7 @@ public:
                   const double& crval1, const double& crval2,
                   const double& crpix1, const double& crpix2,
                   const double& cdelt1, const double& cdelt2);
-    explicit GWcs(const GFitsHDU* hdu);
+    explicit GWcs(const GFitsHDU& hdu);
     GWcs(const GWcs& wcs);
     virtual ~GWcs(void);
@@ -54,8 +54,8 @@ public:
     // Implemented virtual methods
     virtual int         size(void) const;
-    virtual void        read(const GFitsHDU* hdu);
-    virtual void        write(GFitsHDU* hdu) const;
+    virtual void        read(const GFitsHDU& hdu);
+    virtual void        write(GFitsHDU& hdu) const;
     virtual double      omega(const GSkyPixel& pixel) const;
     virtual GSkyDir     pix2dir(const GSkyPixel& pixel) const;
     virtual GSkyPixel   dir2pix(const GSkyDir& dir) const;
@@ -363,11 +363,11 @@ void GEnergies::load(const std::string&amp; filename, const std::string&amp; extname)
     // Open FITS file
     file.open(filename);
-    // Get energies HDU
-    GFitsTable* hdu = file.table(extname);
+    // Get energies table
+    const GFitsTable& table = *file.table(extname);
-    // Read energies from HDU
-    read(hdu);
+    // Read energies from table
+    read(table);
     // Close FITS file
     file.close();
@@ -396,7 +396,7 @@ void GEnergies::save(const std::string&amp; filename, bool clobber,
     GFits file;
     // Write energies to FITS file
-    write(&file, extname);
+    write(file, extname);
     // Save to file
     file.saveto(filename, clobber);
@@ -409,11 +409,11 @@ void GEnergies::save(const std::string&amp; filename, bool clobber,
 /***********************************************************************//**
  * @brief Read energies from FITS table
  *
- * @param[in] hdu Pointer to FITS table.
+ * @param[in] table FITS table.
  *
  * Reads the energies from a FITS table.
  ***************************************************************************/
-void GEnergies::read(const GFitsTable* hdu)
+void GEnergies::read(const GFitsTable& table)
 {
     // Free members
     free_members();
@@ -421,33 +421,28 @@ void GEnergies::read(const GFitsTable* hdu)
     // Initialise attributes
     init_members();
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
+    // Extract number of energy bins in FITS file
+    int num = table.integer("NAXIS2");
-        // Extract number of energy bins in FITS file
-        int num = hdu->integer("NAXIS2");
+    // Continue only if there are energy bins
+    if (num > 0) {
-        // Continue only if there are energy bins
-        if (num > 0) {
-
-            // Get the unit of the energies. If no TUNIT1 header keyword is
-            // found then use MeV
-            std::string unit = "MeV";
-            if (hdu->hascard("TUNIT1")) {
-                unit = hdu->string("TUNIT1");
-            }
+        // Get the unit of the energies. If no TUNIT1 header keyword is
+        // found then use MeV
+        std::string unit = "MeV";
+        if (table.hascard("TUNIT1")) {
+            unit = table.string("TUNIT1");
+        }
-            // Get the column with the name "Energy"
-            const GFitsTableCol* col_energy = (*hdu)["Energy"];
+        // Get the column with the name "Energy"
+        const GFitsTableCol* col_energy = table["Energy"];
-            // Set energies
-            for (int i = 0; i < num; ++i) {
-                append(GEnergy(col_energy->real(i), unit));
-            }
-
-        } // endif: there were energy bins
+        // Set energies
+        for (int i = 0; i < num; ++i) {
+            append(GEnergy(col_energy->real(i), unit));
+        }
-    } // endif: the HDU was valid
+    } // endif: there were energy bins
     // Return
     return;
@@ -457,12 +452,12 @@ void GEnergies::read(const GFitsTable* hdu)
 /***********************************************************************//**
  * @brief Write energies into FITS object
  *
- * @param[in] file Pointer to FITS file.
+ * @param[in] file FITS file.
  * @param[in] extname Energy extension name (defaults to "ENERGIES")
  *
  * Writes energies into FITS object.
  ***************************************************************************/
-void GEnergies::write(GFits* file, const std::string& extname) const
+void GEnergies::write(GFits& file, const std::string& extname) const
 {
     // Set number of energies
     int num = m_energies.size();
@@ -482,7 +477,7 @@ void GEnergies::write(GFits* file, const std::string&amp; extname) const
     table->extname(extname);
     // Write to FITS file
-    file->append(*table);
+    file.append(*table);
     // Free table
     delete table;
@@ -36,7 +36,7 @@
 /* __ Method name definitions ____________________________________________ */
 #define G_CONSTRUCT    "GHealpix::GHealpix(int& ,std::string& ,std::string&)"
-#define G_READ                                    "GHealpix::read(GFitsHDU*)"
+#define G_READ                                    "GHealpix::read(GFitsHDU&)"
 #define G_XY2DIR                               "GHealpix::xy2dir(GSkyPixel&)"
 #define G_DIR2XY2                                "GHealpix::dir2xy(GSkyDir&)"
 #define G_PIX2ANG_RING        "GHealpix::pix2ang_ring(int, double*, double*)"
@@ -133,14 +133,14 @@ GHealpix::GHealpix(const int&amp;         nside,
 /***********************************************************************//**
  * @brief Constructor from FITS HDU table
  *
- * @param[in] hdu Pointer to FITS HDU.
+ * @param[in] hdu FITS HDU.
  ***************************************************************************/
-GHealpix::GHealpix(const GFitsHDU* hdu) : GSkyProjection()
+GHealpix::GHealpix(const GFitsHDU& hdu) : GSkyProjection()
 {
     // Initialise class members
     init_members();
-    // Read Healpix defintion from FITS HDU
+    // Read Healpix defintion from FITS table
     read(hdu);
     // Return
@@ -253,7 +253,7 @@ GHealpix* GHealpix::clone(void) const
 /***********************************************************************//**
  * @brief Read Healpix definiton from FITS header
  *
- * @param[in] hdu FITS HDU containing the Healpix definition.
+ * @param[in] hdu FITS HDU.
  *
  * @exception GException::wcs
  *            Unable to load Healpix definition from HDU.
@@ -262,29 +262,29 @@ GHealpix* GHealpix::clone(void) const
  * @exception GException::wcs_hpx_bad_ordering
  *            Invalid ordering parameter.
  ***************************************************************************/
-void GHealpix::read(const GFitsHDU* hdu)
+void GHealpix::read(const GFitsHDU& hdu)
 {
     // Clear object
     clear();
     // Check if we have a healpix representation
-    if (hdu->string("PIXTYPE") != "HEALPIX") {
+    if (hdu.string("PIXTYPE") != "HEALPIX") {
         throw GException::wcs(G_READ, "HDU does not contain Healpix data");
     }
     // Get pixel ordering
-    std::string order = hdu->string("ORDERING");
+    std::string order = hdu.string("ORDERING");
     // Get coordinate system.
     // First search for HIER_CRD keyword (this has been used in older
     // versions of LAT exposure cubes). If not found then search for standard
     // COORDSYS keyword.
     std::string coords;
-    if (hdu->hascard("HIER_CRD")) {
-        coords = hdu->string("HIER_CRD");
+    if (hdu.hascard("HIER_CRD")) {
+        coords = hdu.string("HIER_CRD");
     }
-    else if (hdu->hascard("COORDSYS")) {
-        coords = hdu->string("COORDSYS");
+    else if (hdu.hascard("COORDSYS")) {
+        coords = hdu.string("COORDSYS");
     }
     // Set coordinate system
@@ -294,7 +294,7 @@ void GHealpix::read(const GFitsHDU* hdu)
     ordering(order);
     // Get Healpix resolution and determine number of pixels and solid angle
-    m_nside      = hdu->integer("NSIDE");
+    m_nside      = hdu.integer("NSIDE");
     m_npface     = m_nside * m_nside;
     m_ncap       = 2 * (m_npface - m_nside);
     m_num_pixels = 12 * m_npface;
@@ -322,25 +322,20 @@ void GHealpix::read(const GFitsHDU* hdu)
  * ORDERING = ordering()
  * COORDSYS = coordsys()
  ***************************************************************************/
-void GHealpix::write(GFitsHDU* hdu) const
+void GHealpix::write(GFitsHDU& hdu) const
 {
-    // Continue only if file pointer is valid
-    if (hdu != NULL) {
-
-        // Set extension name
-        hdu->extname("HEALPIX");
-
-        // Set keywords
-        hdu->card("PIXTYPE",  "HEALPIX",  "HEALPix pixelisation");
-        hdu->card("NSIDE",    nside(),    "HEALPix resolution parameter");
-        hdu->card("FIRSTPIX", 0,          "Index of first pixel");
-        hdu->card("LASTPIX",  npix()-1,   "Index of last pixel");
-        hdu->card("ORDERING", ordering(),
-                  "Pixel ordering scheme, either RING or NESTED");
-        hdu->card("COORDSYS", coordsys(),
-                  "Coordinate system, either EQU or GAL");
-
-    } // endif: HDU was valid
+    // Set extension name
+    hdu.extname("HEALPIX");
+
+    // Set keywords
+    hdu.card("PIXTYPE",  "HEALPIX",  "HEALPix pixelisation");
+    hdu.card("NSIDE",    nside(),    "HEALPix resolution parameter");
+    hdu.card("FIRSTPIX", 0,          "Index of first pixel");
+    hdu.card("LASTPIX",  npix()-1,   "Index of last pixel");
+    hdu.card("ORDERING", ordering(),
+             "Pixel ordering scheme, either RING or NESTED");
+    hdu.card("COORDSYS", coordsys(),
+             "Coordinate system, either EQU or GAL");
     // Return
     return;
@@ -52,7 +52,7 @@
 #define G_DIR2PIX                                "GSkymap::dir2pix(GSkyDir&)"
 #define G_OMEGA1                                       "GSkymap::omega(int&)"
 #define G_OMEGA2                                 "GSkymap::omega(GSkyPixel&)"
-#define G_READ                               "GSkymap::read(const GFitsHDU*)"
+#define G_READ                               "GSkymap::read(const GFitsHDU&)"
 #define G_SET_WCS     "GSkymap::set_wcs(std::string&, std::string&, double&,"\
                               " double&, double&, double&, double&, double&,"\
                                                        " GMatrix&, GVector&)"
@@ -937,21 +937,20 @@ void GSkymap::load(const std::string&amp; filename)
     // Get number of HDUs
     int num = fits.size();
-    // Initialize pointer to HDU and load flag
-    GFitsHDU* hdu    = NULL;
-    bool      loaded = false;
+    // Initialize load flag
+    bool loaded = false;
     // First search for HEALPix extension. We can skip the first extension
     // since this is always an image and a HEALPix map is stored in a
     // binary table
     for (int extno = 1; extno < num; ++extno) {
-        // Get pointer to HDU
-        hdu = fits.at(extno);
+        // Get reference to HDU
+        const GFitsHDU& hdu = *fits.at(extno);
         // If PIXTYPE keyword equals "HEALPIX" then load map
-        if (hdu->hascard("PIXTYPE") && hdu->string("PIXTYPE") == "HEALPIX") {
-            read_healpix(static_cast<const GFitsTable*>(hdu));
+        if (hdu.hascard("PIXTYPE") && hdu.string("PIXTYPE") == "HEALPIX") {
+            read_healpix(static_cast<const GFitsTable&>(hdu));
             loaded = true;
             break;
         }
@@ -963,17 +962,17 @@ void GSkymap::load(const std::string&amp; filename)
     if (!loaded) {
         for (int extno = 0; extno < num; ++extno) {
-            // Get pointer to HDU
-            hdu = fits.at(extno);
+            // Get referene to HDU
+            const GFitsHDU& hdu = *fits.at(extno);
             // Skip if extension is not an image
             if (extno > 0) {
-                if (hdu->string("XTENSION") != "IMAGE")
+                if (hdu.string("XTENSION") != "IMAGE")
                     continue;
             }
             // Load WCS map
-            read_wcs(static_cast<const GFitsImage*>(hdu));
+            read_wcs(static_cast<const GFitsImage&>(hdu));
             loaded = true;
             break;
@@ -1035,39 +1034,32 @@ void GSkymap::save(const std::string&amp; filename, bool clobber) const
  * @brief Read skymap from FITS HDU
  *
  * @param[in] hdu FITS HDU.
- *
- * The method returns an empty skymap of the HDU pointer was not valid.
  ***************************************************************************/
-void GSkymap::read(const GFitsHDU* hdu)
+void GSkymap::read(const GFitsHDU& hdu)
 {
     // Free memory and initialise members
     free_members();
     init_members();
-    // Continue only if HDU pointer is valid
-    if (hdu != NULL) {
+    // Initialize load flag
+    bool loaded = false;
-        // Initialize load flag
-        bool loaded = false;
+    // If PIXTYPE keyword equals "HEALPIX" then load map
+    if (hdu.hascard("PIXTYPE") && hdu.string("PIXTYPE") == "HEALPIX") {
+        read_healpix(static_cast<const GFitsTable&>(hdu));
+        loaded = true;
+    }
-        // If PIXTYPE keyword equals "HEALPIX" then load map
-        if (hdu->hascard("PIXTYPE") && hdu->string("PIXTYPE") == "HEALPIX") {
-            read_healpix(static_cast<const GFitsTable*>(hdu));
+    // ... otherwise try loading as non HEALPix map
+    if (!loaded) {
+
+        // Load only if HDU contains an image
+        if (hdu.exttype() == 0) {
+            read_wcs(static_cast<const GFitsImage&>(hdu));
             loaded = true;
         }
-        // ... otherwise try loading as non HEALPix map
-        if (!loaded) {
-
-            // Load only if HDU contains an image
-            if (hdu->exttype() == 0) {
-                read_wcs(static_cast<const GFitsImage*>(hdu));
-                loaded = true;
-            }
-
-        } // endif
-
-    } // endif: HDU pointer was valid
+    } // endif
     // Return
     return;
@@ -1079,7 +1071,7 @@ void GSkymap::read(const GFitsHDU* hdu)
  *
  * @param[in] file FITS file pointer.
  ***************************************************************************/
-void GSkymap::write(GFits* file) const
+void GSkymap::write(GFits& file) const
 {
     // Continue only if we have data to save
     if (m_proj != NULL) {
@@ -1099,7 +1091,7 @@ void GSkymap::write(GFits* file) const
         // Append HDU to FITS file.
         if (hdu != NULL) {
-            file->append(*hdu);
+            file.append(*hdu);
         }
         // Delete HDU
@@ -1330,7 +1322,7 @@ void GSkymap::set_wcs(const std::string&amp; wcs, const std::string&amp; coords,
 /***********************************************************************//**
  * @brief Read Healpix data from FITS table.
  *
- * @param[in] hdu FITS HDU containing the Healpix data.
+ * @param[in] table FITS table.
  *
  * HEALPix data may be stored in various formats depending on the 
  * application that has writted the data. HEALPix IDL, for example, may
@@ -1339,114 +1331,109 @@ void GSkymap::set_wcs(const std::string&amp; wcs, const std::string&amp; coords,
  * several HEALPix maps into a single column. Alternatively, multiple maps
  * may be stored in multiple columns.
  ***************************************************************************/
-void GSkymap::read_healpix(const GFitsTable* hdu)
+void GSkymap::read_healpix(const GFitsTable& table)
 {
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Determine number of rows and columns in table
-        int nrows = hdu->nrows();
-        int ncols = hdu->ncols();
-        #if defined(G_READ_HEALPIX_DEBUG)
-        std::cout << "nrows=" << nrows << " ncols=" << ncols << std::endl;
-        #endif
-
-        // Allocate Healpix projection
-        m_proj = new GHealpix;
-
-        // Read WCS information from FITS header
-        m_proj->read(hdu);
-
-        // Set number of pixels based on NSIDE parameter
-        m_num_pixels = static_cast<GHealpix*>(m_proj)->npix();
-        #if defined(G_READ_HEALPIX_DEBUG)
-        std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
-        #endif
-
-        // Number of map pixels has to be a multiple of the number of
-        // rows in column
-        if (m_num_pixels % nrows != 0) {
-            throw GException::skymap_bad_size(G_READ_HEALPIX, nrows,
-                                              m_num_pixels);
-        }
+    // Determine number of rows and columns in table
+    int nrows = table.nrows();
+    int ncols = table.ncols();
+    #if defined(G_READ_HEALPIX_DEBUG)
+    std::cout << "nrows=" << nrows << " ncols=" << ncols << std::endl;
+    #endif
-        // Determine vector length for HEALPix data storage
-        int nentry = m_num_pixels / nrows;
-        #if defined(G_READ_HEALPIX_DEBUG)
-        std::cout << "nentry=" << nentry << std::endl;
-        #endif
-
-        // Determine number of maps from the number of maps that fit into
-        // all columns. Only count columns that can fully hold the map.
-        m_num_maps = 0;
-        for (int icol = 0; icol < ncols; ++icol) {
-            const GFitsTableCol* col = (*hdu)[icol];
-            if (col->number() % nentry == 0) {
-                m_num_maps += col->number() / nentry;
-            }
-        }
-        #if defined(G_READ_HEALPIX_DEBUG)
-        std::cout << "m_num_maps=" << m_num_maps << std::endl;
-        #endif
+    // Allocate Healpix projection
+    m_proj = new GHealpix;
-        // Allocate pixels to hold the map
-        alloc_pixels();
+    // Read WCS information from FITS header
+    m_proj->read(table);
+
+    // Set number of pixels based on NSIDE parameter
+    m_num_pixels = static_cast<GHealpix*>(m_proj)->npix();
+    #if defined(G_READ_HEALPIX_DEBUG)
+    std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
+    #endif
-        // Initialise map counter
-        int imap = 0;
+    // Number of map pixels has to be a multiple of the number of
+    // rows in column
+    if (m_num_pixels % nrows != 0) {
+        throw GException::skymap_bad_size(G_READ_HEALPIX, nrows,
+                                          m_num_pixels);
+    }
-        // Loop over all columns
-        for (int icol = 0; icol < ncols; ++icol) {
+    // Determine vector length for HEALPix data storage
+    int nentry = m_num_pixels / nrows;
+    #if defined(G_READ_HEALPIX_DEBUG)
+    std::cout << "nentry=" << nentry << std::endl;
+    #endif
+
+    // Determine number of maps from the number of maps that fit into
+    // all columns. Only count columns that can fully hold the map.
+    m_num_maps = 0;
+    for (int icol = 0; icol < ncols; ++icol) {
+        const GFitsTableCol* col = table[icol];
+        if (col->number() % nentry == 0) {
+            m_num_maps += col->number() / nentry;
+        }
+    }
+    #if defined(G_READ_HEALPIX_DEBUG)
+    std::cout << "m_num_maps=" << m_num_maps << std::endl;
+    #endif
-            // Get next column
-            const GFitsTableCol* col = (*hdu)[icol];
+    // Allocate pixels to hold the map
+    alloc_pixels();
-            // Only consider columns that can fully hold maps
-            if (col->number() % nentry == 0) {
+    // Initialise map counter
+    int imap = 0;
-                // Determine number of maps in column
-                int num = col->number() / nentry;
+    // Loop over all columns
+    for (int icol = 0; icol < ncols; ++icol) {
-                // Loop over all maps in column
-                int inx_start = 0;
-                int inx_end   = nentry;
-                for (int i = 0; i < num; ++i) {
+        // Get next column
+        const GFitsTableCol* col = table[icol];
-                    // Load map
-                    double *ptr = m_pixels + m_num_pixels*imap;
-                    for (int row = 0; row < col->length(); ++row) {
-                        for (int inx = inx_start; inx < inx_end; ++inx) {
-                            *ptr++ = col->real(row,inx);
-                        }
-                    }
-                    #if defined(G_READ_HEALPIX_DEBUG)
-                    std::cout << "Load map=" << imap << " index="
-                              << inx_start << "-" << inx_end << std::endl;
-                    #endif
+        // Only consider columns that can fully hold maps
+        if (col->number() % nentry == 0) {
-                    // Increment index range
-                    inx_start  = inx_end;
-                    inx_end   += nentry;
+            // Determine number of maps in column
+            int num = col->number() / nentry;
-                    // Increment map counter
-                    imap++;
+            // Loop over all maps in column
+            int inx_start = 0;
+            int inx_end   = nentry;
+            for (int i = 0; i < num; ++i) {
-                    // Break if we have loaded all maps
-                    if (imap >= m_num_maps) {
-                        break;
+                // Load map
+                double *ptr = m_pixels + m_num_pixels*imap;
+                for (int row = 0; row < col->length(); ++row) {
+                    for (int inx = inx_start; inx < inx_end; ++inx) {
+                        *ptr++ = col->real(row,inx);
                     }
+                }
+                #if defined(G_READ_HEALPIX_DEBUG)
+                std::cout << "Load map=" << imap << " index="
+                          << inx_start << "-" << inx_end << std::endl;
+                #endif
-                } // endfor: looped over all maps in column
-            } // endif: column could fully hold maps
+                // Increment index range
+                inx_start  = inx_end;
+                inx_end   += nentry;
-            // Break if we have loaded all maps
-            if (imap >= m_num_maps) {
-                break;
-            }
+                // Increment map counter
+                imap++;
-        } // endfor: looped over all columns
+                // Break if we have loaded all maps
+                if (imap >= m_num_maps) {
+                    break;
+                }
+
+            } // endfor: looped over all maps in column
+        } // endif: column could fully hold maps
+
+        // Break if we have loaded all maps
+        if (imap >= m_num_maps) {
+            break;
+        }
-    } // endif: HDU was valid
+    } // endfor: looped over all columns
     // Return
     return;
@@ -1456,72 +1443,67 @@ void GSkymap::read_healpix(const GFitsTable* hdu)
 /***********************************************************************//**
  * @brief Read WCS image from FITS HDU
  *
- * @param[in] hdu FITS HDU containing the WCS image.
+ * @param[in] image FITS image.
  *
  * @exception GException::skymap_bad_image_dim
  *            WCS image has invalid dimension (naxis=2 or 3).
  ***************************************************************************/
-void GSkymap::read_wcs(const GFitsImage* hdu)
+void GSkymap::read_wcs(const GFitsImage& image)
 {
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
+    // Allocate WCS
+    alloc_wcs(image);
-        // Allocate WCS
-        alloc_wcs(hdu);
+    // Read projection information from FITS header
+    m_proj->read(image);
-        // Read projection information from FITS header
-        m_proj->read(hdu);
+    // Extract map dimension and number of maps from image
+    if (image.naxis() == 2) {
+        m_num_x    = image.naxes(0);
+        m_num_y    = image.naxes(1);
+        m_num_maps = 1;
+    }
+    else if (image.naxis() >= 3) {
+        m_num_x    = image.naxes(0);
+        m_num_y    = image.naxes(1);
+        m_num_maps = image.naxes(2);
+    }
+    else {
+        throw GException::skymap_bad_image_dim(G_READ_WCS, image.naxis());
+    }
+    #if defined(G_READ_WCS_DEBUG)
+    std::cout << "m_num_x=" << m_num_x << std::endl;
+    std::cout << "m_num_y=" << m_num_y << std::endl;
+    std::cout << "m_num_maps=" << m_num_maps << std::endl;
+    #endif
-        // Extract map dimension and number of maps from image
-        if (hdu->naxis() == 2) {
-            m_num_x    = hdu->naxes(0);
-            m_num_y    = hdu->naxes(1);
-            m_num_maps = 1;
-        }
-        else if (hdu->naxis() >= 3) {
-            m_num_x    = hdu->naxes(0);
-            m_num_y    = hdu->naxes(1);
-            m_num_maps = hdu->naxes(2);
-        }
-        else {
-            throw GException::skymap_bad_image_dim(G_READ_WCS, hdu->naxis());
-        }
-        #if defined(G_READ_WCS_DEBUG)
-        std::cout << "m_num_x=" << m_num_x << std::endl;
-        std::cout << "m_num_y=" << m_num_y << std::endl;
-        std::cout << "m_num_maps=" << m_num_maps << std::endl;
-        #endif
-
-        // Compute number of pixels
-        m_num_pixels = m_num_x * m_num_y;
-        #if defined(G_READ_WCS_DEBUG)
-        std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
-        #endif
-
-        // Allocate pixels to hold the map
-        alloc_pixels();
+    // Compute number of pixels
+    m_num_pixels = m_num_x * m_num_y;
+    #if defined(G_READ_WCS_DEBUG)
+    std::cout << "m_num_pixels=" << m_num_pixels << std::endl;
+    #endif
-        // Read image
-        if (hdu->naxis() == 2) {
-            double* ptr = m_pixels;
-            for (int iy = 0; iy < m_num_y; ++iy) {
-                for (int ix = 0; ix < m_num_x; ++ix) {
-                    *ptr++ = hdu->pixel(ix,iy);
-                }
+    // Allocate pixels to hold the map
+    alloc_pixels();
+
+    // Read image
+    if (image.naxis() == 2) {
+        double* ptr = m_pixels;
+        for (int iy = 0; iy < m_num_y; ++iy) {
+            for (int ix = 0; ix < m_num_x; ++ix) {
+                *ptr++ = image.pixel(ix,iy);
             }
         }
-        else {
-            double* ptr = m_pixels;
-            for (int imap = 0; imap < m_num_maps; ++imap) {
-                for (int iy = 0; iy < m_num_y; ++iy) {
-                    for (int ix = 0; ix < m_num_x; ++ix) {
-                        *ptr++ = hdu->pixel(ix,iy,imap);
-                    }
+    }
+    else {
+        double* ptr = m_pixels;
+        for (int imap = 0; imap < m_num_maps; ++imap) {
+            for (int iy = 0; iy < m_num_y; ++iy) {
+                for (int ix = 0; ix < m_num_x; ++ix) {
+                    *ptr++ = image.pixel(ix,iy,imap);
                 }
             }
         }
-
-    } // endif: HDU was valid
+    }
     // Return
     return;
@@ -1531,7 +1513,7 @@ void GSkymap::read_wcs(const GFitsImage* hdu)
 /***********************************************************************//**
  * @brief Allocate WCS class
  *
- * @param[in] hdu FITS HDU containing the WCS image.
+ * @param[in] image FITS image.
  *
  * @exception GException::fits_key_not_found
  *            Unable to find required FITS header keyword.
@@ -1540,39 +1522,34 @@ void GSkymap::read_wcs(const GFitsImage* hdu)
  * @exception GException::wcs_invalid
  *            WCS projection of FITS file not supported by GammaLib.
  ***************************************************************************/
-void GSkymap::alloc_wcs(const GFitsImage* hdu)
+void GSkymap::alloc_wcs(const GFitsImage& image)
 {
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Get standard keywords
-        std::string ctype1 = hdu->string("CTYPE1");
-        std::string ctype2 = hdu->string("CTYPE2");
-
-        // Extract projection type
-        std::string xproj = ctype1.substr(5,3);
-        std::string yproj = ctype2.substr(5,3);
-
-        // Check that projection type is identical on both axes
-        if (xproj != yproj) {
-            throw GException::skymap_bad_ctype(G_ALLOC_WCS,
-                                               ctype1, ctype2);
-        }
-
-        // Allocate WCS registry
-        GWcsRegistry registry;
-        
-        // Allocate projection from registry
-        m_proj = registry.alloc(xproj);
-        
-        // Signal if projection type is not known
-        if (m_proj == NULL) {
-            std::string message = "Projection code not known. "
-                                  "Should be one of "+registry.list()+".";
-            throw GException::wcs_invalid(G_ALLOC_WCS, xproj, message);
-        }
+    // Get standard keywords
+    std::string ctype1 = image.string("CTYPE1");
+    std::string ctype2 = image.string("CTYPE2");
+
+    // Extract projection type
+    std::string xproj = ctype1.substr(5,3);
+    std::string yproj = ctype2.substr(5,3);
+
+    // Check that projection type is identical on both axes
+    if (xproj != yproj) {
+        throw GException::skymap_bad_ctype(G_ALLOC_WCS,
+                                           ctype1, ctype2);
+    }
-    } // endif: HDU was valid
+    // Allocate WCS registry
+    GWcsRegistry registry;
+    
+    // Allocate projection from registry
+    m_proj = registry.alloc(xproj);
+    
+    // Signal if projection type is not known
+    if (m_proj == NULL) {
+        std::string message = "Projection code not known. "
+                              "Should be one of "+registry.list()+".";
+        throw GException::wcs_invalid(G_ALLOC_WCS, xproj, message);
+    }
     // Return
     return;
@@ -1626,7 +1603,7 @@ GFitsBinTable* GSkymap::create_healpix_hdu(void) const
     hdu->extname("HEALPIX");
     // If we have WCS information then write into FITS header
-    if (m_proj != NULL) m_proj->write(hdu);
+    if (m_proj != NULL) m_proj->write(*hdu);
     // Set additional keywords
     hdu->card("NBRBINS", m_num_maps, "Number of HEALPix maps");
@@ -1693,7 +1670,7 @@ GFitsImageDouble* GSkymap::create_wcs_hdu(void) const
     hdu->extname("IMAGE");
     // If we have sky projection information then write into FITS header
-    if (m_proj != NULL) m_proj->write(hdu);
+    if (m_proj != NULL) m_proj->write(*hdu);
     // Set additional keywords
     //TODO
@@ -36,7 +36,7 @@
 #include "GWcs.hpp"
 /* __ Method name definitions ____________________________________________ */
-#define G_READ                                        "GWcs::read(GFitsHDU*)"
+#define G_READ                                        "GWcs::read(GFitsHDU&)"
 #define G_OMEGA                                           "GWcs::omega(int&)"
 #define G_PIX2DIR                                       "GWcs::pix2dir(int&)"
 #define G_DIR2PIX                                   "GWcs::dir2pix(GSkyDir&)"
@@ -121,7 +121,7 @@ GWcs::GWcs(const std::string&amp; coords,
  *
  * @param[in] hdu FITS HDU.
  ***************************************************************************/
-GWcs::GWcs(const GFitsHDU* hdu) : GSkyProjection()
+GWcs::GWcs(const GFitsHDU& hdu) : GSkyProjection()
 {
     // Initialise class members
     init_members();
@@ -210,93 +210,88 @@ GWcs&amp; GWcs::operator=(const GWcs&amp; wcs)
 /***********************************************************************//**
  * @brief Read WCS definiton from FITS header.
  *
- * @param[in] hdu FITS HDU containing the WCS definition.
+ * @param[in] hdu FITS HDU.
  *
  * @exception GException::wcs_bad_coords
  *            Coordinate system is not of valid type.
  *
  * This method reads the WCS definition from the FITS header.
  ***************************************************************************/
-void GWcs::read(const GFitsHDU* hdu)
+void GWcs::read(const GFitsHDU& hdu)
 {
     // Clear object
     clear();
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Get standard keywords
-        std::string ctype1 = hdu->string("CTYPE1");
-        std::string ctype2 = hdu->string("CTYPE2");
-        double      crval1 = hdu->real("CRVAL1");
-        double      crval2 = hdu->real("CRVAL2");
-        double      crpix1 = hdu->real("CRPIX1");
-        double      crpix2 = hdu->real("CRPIX2");
-        double      cdelt1 = hdu->real("CDELT1");
-        double      cdelt2 = hdu->real("CDELT2");
-
-        // Determine coordinate system
-        std::string coords;
-        std::string xcoord = ctype1.substr(0,4);
-        std::string ycoord = ctype2.substr(0,4);
-        if (xcoord == "RA--" && ycoord == "DEC-") {
-            m_lng  = 0;
-            m_lat  = 1;
-            coords = "EQU";
-        }
-        else if (xcoord == "DEC-" && ycoord == "RA--") {
-            m_lng  = 1;
-            m_lat  = 0;
-            coords = "EQU";
-        }
-        else if (xcoord == "GLON" && ycoord == "GLAT") {
-            m_lng  = 0;
-            m_lat  = 1;
-            coords = "GAL";
-        }
-        else if (xcoord == "GLAT" && ycoord == "GLON") {
-            m_lng  = 1;
-            m_lat  = 0;
-            coords = "GAL";
-        }
-        else if (xcoord == "ELON" && ycoord == "ELAT") {
-            m_lng  = 0;
-            m_lat  = 1;
-            coords = "ECL";
-        }
-        else if (xcoord == "ELAT" && ycoord == "ELON") {
-            m_lng  = 1;
-            m_lat  = 0;
-            coords = "ECL";
-        }
-        else if (xcoord == "HLON" && ycoord == "HLAT") {
-            m_lng  = 0;
-            m_lat  = 1;
-            coords = "HEL";
-        }
-        else if (xcoord == "HLAT" && ycoord == "HLON") {
-            m_lng  = 1;
-            m_lat  = 0;
-            coords = "HEL";
-        }
-        else if (xcoord == "SLON" && ycoord == "SLAT") {
-            m_lng  = 0;
-            m_lat  = 1;
-            coords = "SGL";
-        }
-        else if (xcoord == "SLAT" && ycoord == "SLON") {
-            m_lng  = 1;
-            m_lat  = 0;
-            coords = "SGL";
-        }
-        else {
-            throw GException::wcs_bad_coords(G_READ, coordsys());
-        }
-
-        // Set standard parameters
-        set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
+    // Get standard keywords
+    std::string ctype1 = hdu.string("CTYPE1");
+    std::string ctype2 = hdu.string("CTYPE2");
+    double      crval1 = hdu.real("CRVAL1");
+    double      crval2 = hdu.real("CRVAL2");
+    double      crpix1 = hdu.real("CRPIX1");
+    double      crpix2 = hdu.real("CRPIX2");
+    double      cdelt1 = hdu.real("CDELT1");
+    double      cdelt2 = hdu.real("CDELT2");
+
+    // Determine coordinate system
+    std::string coords;
+    std::string xcoord = ctype1.substr(0,4);
+    std::string ycoord = ctype2.substr(0,4);
+    if (xcoord == "RA--" && ycoord == "DEC-") {
+        m_lng  = 0;
+        m_lat  = 1;
+        coords = "EQU";
+    }
+    else if (xcoord == "DEC-" && ycoord == "RA--") {
+        m_lng  = 1;
+        m_lat  = 0;
+        coords = "EQU";
+    }
+    else if (xcoord == "GLON" && ycoord == "GLAT") {
+        m_lng  = 0;
+        m_lat  = 1;
+        coords = "GAL";
+    }
+    else if (xcoord == "GLAT" && ycoord == "GLON") {
+        m_lng  = 1;
+        m_lat  = 0;
+        coords = "GAL";
+    }
+    else if (xcoord == "ELON" && ycoord == "ELAT") {
+        m_lng  = 0;
+        m_lat  = 1;
+        coords = "ECL";
+    }
+    else if (xcoord == "ELAT" && ycoord == "ELON") {
+        m_lng  = 1;
+        m_lat  = 0;
+        coords = "ECL";
+    }
+    else if (xcoord == "HLON" && ycoord == "HLAT") {
+        m_lng  = 0;
+        m_lat  = 1;
+        coords = "HEL";
+    }
+    else if (xcoord == "HLAT" && ycoord == "HLON") {
+        m_lng  = 1;
+        m_lat  = 0;
+        coords = "HEL";
+    }
+    else if (xcoord == "SLON" && ycoord == "SLAT") {
+        m_lng  = 0;
+        m_lat  = 1;
+        coords = "SGL";
+    }
+    else if (xcoord == "SLAT" && ycoord == "SLON") {
+        m_lng  = 1;
+        m_lat  = 0;
+        coords = "SGL";
+    }
+    else {
+        throw GException::wcs_bad_coords(G_READ, coordsys());
+    }
-    } // endif: HDU was valid
+    // Set standard parameters
+    set_members(coords, crval1, crval2, crpix1, crpix2, cdelt1, cdelt2);
     // Return
     return;
@@ -311,122 +306,115 @@ void GWcs::read(const GFitsHDU* hdu)
  * This method writes the World Coordinate System definition into the FITS
  * HDU header.
  *
- * The method does nothing if the HDU is not valid.
- *
  * This method has been adapted from wcshdr.c::wcshdo().
  ***************************************************************************/
-void GWcs::write(GFitsHDU* hdu) const
+void GWcs::write(GFitsHDU& hdu) const
 {
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Write reference pixel coordinates
-        for (int i = 0; i < m_naxis; ++i) {
-            std::string keyname = "CRPIX"+gammalib::str(i+1);
-            std::string comment = "Pixel coordinate of reference point (starting from 1)";
-            hdu->card(keyname, m_crpix.at(i), comment);
-        }
-        
-        //TODO: Write linear transformation matrix
-        
-        // Write coordinate increment at reference point
-        for (int i = 0; i < m_naxis; ++i) {
-            std::string keyname = "CDELT"+gammalib::str(i+1);
-            std::string comment;
-            if (m_cunit.at(i).length() > 0) {
-                comment += "["+gammalib::strip_whitespace(m_cunit.at(i))+"] ";
-            }
-            comment += "Coordinate increment at reference point";
-            hdu->card(keyname, m_cdelt.at(i), comment);
-        }
-        
-        // Write units of coordinate increment and reference value
-        for (int i = 0; i < m_naxis; ++i) {
-            if (m_cunit.at(i).length() > 0) {
-                std::string keyname = "CUNIT"+gammalib::str(i+1);
-                std::string comment = "Units of coordinate increment and value";
-                hdu->card(keyname, m_cunit.at(0), comment);
-            }
-        }
-        
-        // Write coordinate type
-        wcs_set_ctype();
-        for (int i = 0; i < m_naxis; ++i) {
-            if (i == m_lng) {
-                std::string keyname = "CTYPE"+gammalib::str(i+1);
-                hdu->card(keyname, m_ctype.at(i), m_ctype_c.at(i));
-            }
-            if (i == m_lat) {
-                std::string keyname = "CTYPE"+gammalib::str(i+1);
-                hdu->card(keyname, m_ctype.at(i), m_ctype_c.at(i));
-            }
-            if (i == m_spec) {
-                std::string keyname = "CTYPE"+gammalib::str(i+1);
-                hdu->card(keyname, m_ctype.at(i), m_ctype_c.at(i));
-            }
-        }
-
-        // Write coordinate value at reference point
-        for (int i = 0; i < m_naxis; ++i) {
-            std::string keyname = "CRVAL"+gammalib::str(i+1);
-            std::string comment;
-            if (m_cunit.at(i).length() > 0) {
-                comment += "["+gammalib::strip_whitespace(m_cunit.at(i))+"] ";
-            }
-            comment += "Coordinate value at reference point";
-            hdu->card(keyname, m_crval.at(i), comment);
-        }
-        
-        //TODO: Parameter values
-        hdu->card("CROTA2", 0.0, "[deg] Rotation Angle"); // Old style, use PV instead
-        //hdu->card("PV2_1",   0.0,          "Projection parameter 1");
-        //hdu->card("PV2_2",   0.0,          "Projection parameter 2");
-        
-        // Celestial and spectral transformation parameters
-        if (!undefined(m_lonpole)) {
-            hdu->card("LONPOLE", m_lonpole, "[deg] Native longitude of celestial pole");
-        }
-        if (!undefined(m_latpole)) {
-            hdu->card("LATPOLE", m_latpole, "[deg] Native latitude of celestial pole");
+    // Write reference pixel coordinates
+    for (int i = 0; i < m_naxis; ++i) {
+        std::string keyname = "CRPIX"+gammalib::str(i+1);
+        std::string comment = "Pixel coordinate of reference point (starting from 1)";
+        hdu.card(keyname, m_crpix.at(i), comment);
+    }
+    
+    //TODO: Write linear transformation matrix
+    
+    // Write coordinate increment at reference point
+    for (int i = 0; i < m_naxis; ++i) {
+        std::string keyname = "CDELT"+gammalib::str(i+1);
+        std::string comment;
+        if (m_cunit.at(i).length() > 0) {
+            comment += "["+gammalib::strip_whitespace(m_cunit.at(i))+"] ";
         }
-        if (!undefined(m_restfrq)) {
-            hdu->card("RESTFRQ", m_restfrq, "[Hz] Line rest frequency");
+        comment += "Coordinate increment at reference point";
+        hdu.card(keyname, m_cdelt.at(i), comment);
+    }
+    
+    // Write units of coordinate increment and reference value
+    for (int i = 0; i < m_naxis; ++i) {
+        if (m_cunit.at(i).length() > 0) {
+            std::string keyname = "CUNIT"+gammalib::str(i+1);
+            std::string comment = "Units of coordinate increment and value";
+            hdu.card(keyname, m_cunit.at(0), comment);
         }
-        if (!undefined(m_restwav)) {
-            hdu->card("RESTWAV", m_restwav, "[Hz] Line rest wavelength");
+    }
+    
+    // Write coordinate type
+    wcs_set_ctype();
+    for (int i = 0; i < m_naxis; ++i) {
+        if (i == m_lng) {
+            std::string keyname = "CTYPE"+gammalib::str(i+1);
+            hdu.card(keyname, m_ctype.at(i), m_ctype_c.at(i));
         }
-        
-        // Equatorial coordinate system type
-        if (m_radesys.length() > 0) {
-            hdu->card("RADESYS", m_radesys, "Equatorial coordinate system");
+        if (i == m_lat) {
+            std::string keyname = "CTYPE"+gammalib::str(i+1);
+            hdu.card(keyname, m_ctype.at(i), m_ctype_c.at(i));
         }
-        
-        // Equinox of equatorial coordinate system
-        if (!undefined(m_equinox)) {
-            hdu->card("EQUINOX", m_equinox, "[yr] Equinox of equatorial coordinates");
+        if (i == m_spec) {
+            std::string keyname = "CTYPE"+gammalib::str(i+1);
+            hdu.card(keyname, m_ctype.at(i), m_ctype_c.at(i));
         }
-        
-        //TODO: Reference frame of spectral coordinates
-        
-        //TODO: Reference frame of spectral observation
-        
-        //TODO: Observer's velocity towards source
-        
-        //TODO: Reference frame of source redshift
-        
-        //TODO: Redshift of the source
-        
-        //TODO: Observatory coordinates
-        
-        //TODO: MJD of observation
-        
-        //TODO: MJD mid-observation time
-        
-        //TODO: ISO-8601 date corresponding to MJD-OBS
-        
-        //TODO: ISO-8601 date corresponding to MJD-AVG
+    }
-    } // endif: HDU was valid
+    // Write coordinate value at reference point
+    for (int i = 0; i < m_naxis; ++i) {
+        std::string keyname = "CRVAL"+gammalib::str(i+1);
+        std::string comment;
+        if (m_cunit.at(i).length() > 0) {
+            comment += "["+gammalib::strip_whitespace(m_cunit.at(i))+"] ";
+        }
+        comment += "Coordinate value at reference point";
+        hdu.card(keyname, m_crval.at(i), comment);
+    }
+    
+    //TODO: Parameter values
+    hdu.card("CROTA2", 0.0, "[deg] Rotation Angle"); // Old style, use PV instead
+    //hdu.card("PV2_1",   0.0,          "Projection parameter 1");
+    //hdu.card("PV2_2",   0.0,          "Projection parameter 2");
+    
+    // Celestial and spectral transformation parameters
+    if (!undefined(m_lonpole)) {
+        hdu.card("LONPOLE", m_lonpole, "[deg] Native longitude of celestial pole");
+    }
+    if (!undefined(m_latpole)) {
+        hdu.card("LATPOLE", m_latpole, "[deg] Native latitude of celestial pole");
+    }
+    if (!undefined(m_restfrq)) {
+        hdu.card("RESTFRQ", m_restfrq, "[Hz] Line rest frequency");
+    }
+    if (!undefined(m_restwav)) {
+        hdu.card("RESTWAV", m_restwav, "[Hz] Line rest wavelength");
+    }
+    
+    // Equatorial coordinate system type
+    if (m_radesys.length() > 0) {
+        hdu.card("RADESYS", m_radesys, "Equatorial coordinate system");
+    }
+    
+    // Equinox of equatorial coordinate system
+    if (!undefined(m_equinox)) {
+        hdu.card("EQUINOX", m_equinox, "[yr] Equinox of equatorial coordinates");
+    }
+    
+    //TODO: Reference frame of spectral coordinates
+    
+    //TODO: Reference frame of spectral observation
+    
+    //TODO: Observer's velocity towards source
+    
+    //TODO: Reference frame of source redshift
+    
+    //TODO: Redshift of the source
+    
+    //TODO: Observatory coordinates
+    
+    //TODO: MJD of observation
+    
+    //TODO: MJD mid-observation time
+    
+    //TODO: ISO-8601 date corresponding to MJD-OBS
+    
+    //TODO: ISO-8601 date corresponding to MJD-AVG
     // Return
     return;
@@ -259,7 +259,7 @@ void GArf::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 data
     read(table);
@@ -303,44 +303,39 @@ void GArf::save(const std::string&amp; filename, const bool&amp; clobber) const
 /***********************************************************************//**
  * @brief Read Auxiliary Response File
  *
- * @param[in] hdu ARF FITS table.
+ * @param[in] table ARF FITS table.
  ***************************************************************************/
-void GArf::read(const GFitsTable* hdu)
+void GArf::read(const GFitsTable& table)
 {
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Get pointer to data columns
-        const GFitsTableCol* energy_lo = (*hdu)["ENERG_LO"];
-        const GFitsTableCol* energy_hi = (*hdu)["ENERG_HI"];
-        const GFitsTableCol* specresp  = (*hdu)["SPECRESP"];
-
-        // Determine effective area conversion factor. Internal
-        // units are cm^2
-        std::string u_specresp = gammalib::tolower(gammalib::strip_whitespace(specresp->unit()));
-        double      c_specresp = 1.0;
-        if (u_specresp == "m^2" || u_specresp == "m2") {
-            c_specresp = 10000.0;
-        }
+    // Get pointer to data columns
+    const GFitsTableCol* energy_lo = table["ENERG_LO"];
+    const GFitsTableCol* energy_hi = table["ENERG_HI"];
+    const GFitsTableCol* specresp  = table["SPECRESP"];
+
+    // Determine effective area conversion factor. Internal
+    // units are cm^2
+    std::string u_specresp = gammalib::tolower(gammalib::strip_whitespace(specresp->unit()));
+    double      c_specresp = 1.0;
+    if (u_specresp == "m^2" || u_specresp == "m2") {
+        c_specresp = 10000.0;
+    }
-        // Extract number of energy bins
-        int num = energy_lo->length();
+    // Extract number of energy bins
+    int num = energy_lo->length();
-        // Set energy bins
-        for (int i = 0; i < num; ++i) {
+    // Set energy bins
+    for (int i = 0; i < num; ++i) {
-            // Append energy bin
-            GEnergy emin(energy_lo->real(i), energy_lo->unit());
-            GEnergy emax(energy_hi->real(i), energy_hi->unit());
-            m_ebounds.append(emin, emax);
-
-            // Append effective area value
-            double aeff = specresp->real(i) * c_specresp;
-            m_specresp.push_back(aeff);
+        // Append energy bin
+        GEnergy emin(energy_lo->real(i), energy_lo->unit());
+        GEnergy emax(energy_hi->real(i), energy_hi->unit());
+        m_ebounds.append(emin, emax);
-        } // endfor: looped over energy bins
+        // Append effective area value
+        double aeff = specresp->real(i) * c_specresp;
+        m_specresp.push_back(aeff);
-    } // endif: HDU was valid
+    } // endfor: looped over energy bins
     // Return
     return;
@@ -339,7 +339,7 @@ void GPha::load(const std::string&amp; filename)
     GFits fits(filename);
     // Get PHA table
-    GFitsTable* pha = fits.table("SPECTRUM");
+    const GFitsTable& pha = *fits.table("SPECTRUM");
     // Read PHA data
     read(pha);
@@ -394,44 +394,39 @@ void GPha::save(const std::string&amp; filename, const bool&amp; clobber) const
 /***********************************************************************//**
  * @brief Read Pulse Height Analyzer spectrum
  *
- * @param[in] hdu PHA FITS table.
+ * @param[in] table PHA FITS table.
  *
  * @exception GException::invalid_value
  *            Mismatch between PHA file and energy boundaries.
  ***************************************************************************/
-void GPha::read(const GFitsTable* hdu)
+void GPha::read(const GFitsTable& table)
 {
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Get data column
-        const GFitsTableCol* col_data = (*hdu)["COUNTS"];
-
-        // Extract number of channels in FITS file
-        int length = hdu->integer("NAXIS2");
-
-        // Check whether column length is okay
-        if (m_ebounds.size() > 0) {
-            if (m_ebounds.size() != length) {
-                std::string msg = "Number of channels in PHA file ("
-                                  ""+gammalib::str(length)+") mismatches the"
-                                  " number of energy boundaries ("
-                                  ""+gammalib::str(m_ebounds.size())+") that"
-                                  " are defined in the GPha instance.\n"
-                                  "Please define the correct energy boundaries.";
-                throw GException::invalid_value(G_READ, msg);
-            }
+    // Get data column
+    const GFitsTableCol* col_data = table["COUNTS"];
+
+    // Extract number of channels in FITS file
+    int length = table.integer("NAXIS2");
+
+    // Check whether column length is okay
+    if (m_ebounds.size() > 0) {
+        if (m_ebounds.size() != length) {
+            std::string msg = "Number of channels in PHA file ("
+                              ""+gammalib::str(length)+") mismatches the"
+                              " number of energy boundaries ("
+                              ""+gammalib::str(m_ebounds.size())+") that"
+                              " are defined in the GPha instance.\n"
+                              "Please define the correct energy boundaries.";
+            throw GException::invalid_value(G_READ, msg);
         }
+    }
-        // Initialize spectrum
-        m_counts.assign(length, 0.0);
-
-        // Copy data
-        for (int i = 0; i < length; ++i) {
-            m_counts[i] = col_data->real(i);
-        }
+    // Initialize spectrum
+    m_counts.assign(length, 0.0);
-    } // endif: HDU was valid
+    // Copy data
+    for (int i = 0; i < length; ++i) {
+        m_counts[i] = col_data->real(i);
+    }
     // Return
     return;
@@ -276,7 +276,7 @@ void GRmf::load(const std::string&amp; filename)
     m_ebds_measured.load(filename);
     // Get RMF table
-    GFitsTable* table = file.table("MATRIX");
+    const GFitsTable& table = *file.table("MATRIX");
     // Read RMF data
     read(table);
@@ -320,58 +320,53 @@ void GRmf::save(const std::string&amp; filename, const bool&amp; clobber) const
 /***********************************************************************//**
  * @brief Read Redistribution Matrix File
  *
- * @param[in] hdu RMF FITS table.
+ * @param[in] table RMF FITS table.
  ***************************************************************************/
-void GRmf::read(const GFitsTable* hdu)
+void GRmf::read(const GFitsTable& table)
 {
-    // Continue only if HDU is valid
-    if (hdu != NULL) {
-
-        // Initialise members
-        m_ebds_true.clear();
+    // Initialise members
+    m_ebds_true.clear();
-        // Get pointer to data columns
-        const GFitsTableCol* energy_lo = (*hdu)["ENERG_LO"];
-        const GFitsTableCol* energy_hi = (*hdu)["ENERG_HI"];
-        const GFitsTableCol* n_grp     = (*hdu)["N_GRP"];
-        const GFitsTableCol* f_chan    = (*hdu)["F_CHAN"];
-        const GFitsTableCol* n_chan    = (*hdu)["N_CHAN"];
-        const GFitsTableCol* matrix    = (*hdu)["MATRIX"];
+    // Get pointer to data columns
+    const GFitsTableCol* energy_lo = table["ENERG_LO"];
+    const GFitsTableCol* energy_hi = table["ENERG_HI"];
+    const GFitsTableCol* n_grp     = table["N_GRP"];
+    const GFitsTableCol* f_chan    = table["F_CHAN"];
+    const GFitsTableCol* n_chan    = table["N_CHAN"];
+    const GFitsTableCol* matrix    = table["MATRIX"];
-        // Set matrix rows and columns
-        int rows    = energy_lo->length();
-        int columns = m_ebds_measured.size();
+    // Set matrix rows and columns
+    int rows    = energy_lo->length();
+    int columns = m_ebds_measured.size();
-        // Initialize matrix
-        m_matrix = GMatrixSparse(rows, columns);
+    // Initialize matrix
+    m_matrix = GMatrixSparse(rows, columns);
-        // Set true energy bins
-        for (int itrue = 0; itrue < rows; ++itrue) {
+    // Set true energy bins
+    for (int itrue = 0; itrue < rows; ++itrue) {
-            // Append energy bin
-            GEnergy emin(energy_lo->real(itrue), energy_lo->unit());
-            GEnergy emax(energy_hi->real(itrue), energy_hi->unit());
-            m_ebds_true.append(emin, emax);
-
-            // Loop over groups
-            int icolumn = 0;
-            int ngroups = n_grp->integer(itrue);
-            for (int igroup = 0; igroup < ngroups; ++igroup) {
-
-                // Get start column index and number of columns
-                int imeasured = f_chan->integer(itrue, igroup);
-                int nvalues   = n_chan->integer(itrue, igroup);
-
-                // Get values
-                for (int i = 0; i < nvalues; ++i, ++imeasured, ++icolumn) {
-                    m_matrix(itrue, imeasured) = matrix->real(itrue, icolumn);
-                }
-
-            } // endfor: looped over groups
+        // Append energy bin
+        GEnergy emin(energy_lo->real(itrue), energy_lo->unit());
+        GEnergy emax(energy_hi->real(itrue), energy_hi->unit());
+        m_ebds_true.append(emin, emax);
+
+        // Loop over groups
+        int icolumn = 0;
+        int ngroups = n_grp->integer(itrue);
+        for (int igroup = 0; igroup < ngroups; ++igroup) {
+
+            // Get start column index and number of columns
+            int imeasured = f_chan->integer(itrue, igroup);
+            int nvalues   = n_chan->integer(itrue, igroup);
+
+            // Get values
+            for (int i = 0; i < nvalues; ++i, ++imeasured, ++icolumn) {
+                m_matrix(itrue, imeasured) = matrix->real(itrue, icolumn);
+            }
-        } // endfor: looped over true energy bins
+        } // endfor: looped over groups
-    } // endif: HDU was valid
+    } // endfor: looped over true energy bins
     // Return
     return;