From e821ba0f99688d28d559ce98df3fa4818aafdfe2 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sat, 7 Dec 2013 00:28:22 +0100 Subject: [PATCH] The FITS read() and write() methods now take reference instead of pointers. --- NEWS | 16 +++++++++++++--- include/GArf.hpp | 5 +++-- include/GEnergies.hpp | 4 ++-- include/GHealpix.hpp | 6 +++--- include/GPha.hpp | 5 +++-- include/GPhotons.hpp | 28 ++++++++++++++++++++++++++-- include/GRmf.hpp | 5 +++-- include/GSkyProjection.hpp | 4 ++-- include/GSkymap.hpp | 10 +++++----- include/GTimes.hpp | 28 ++++++++++++++++++++++++++-- include/GWcs.hpp | 6 +++--- inst/com/src/GCOMEventCube.cpp | 22 +++++++++++----------- inst/com/src/GCOMObservation.cpp | 16 ++++++++-------- inst/cta/pyext/GCTAInstDir.i | 5 ++--- inst/cta/src/GCTAEventCube.cpp | 4 ++-- inst/lat/src/GLATEventCube.cpp | 4 ++-- inst/lat/src/GLATLtCubeMap.cpp | 2 +- pyext/GArf.i | 5 +++-- pyext/GEnergies.i | 4 ++-- pyext/GHealpix.i | 6 +++--- pyext/GPha.i | 5 +++-- pyext/GRmf.i | 5 +++-- pyext/GSkyProjection.i | 4 ++-- pyext/GSkymap.i | 4 ++-- pyext/GWcs.i | 6 +++--- src/obs/GEnergies.cpp | 59 +++++++++++++++++++++++++++-------------------------------- src/sky/GHealpix.cpp | 57 ++++++++++++++++++++++++++------------------------------- src/sky/GSkymap.cpp | 403 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- src/sky/GWcs.cpp | 358 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- src/xspec/GArf.cpp | 59 +++++++++++++++++++++++++++-------------------------------- src/xspec/GPha.cpp | 57 ++++++++++++++++++++++++++------------------------------- src/xspec/GRmf.cpp | 83 +++++++++++++++++++++++++++++++++++++++-------------------------------------------- 32 files changed, 644 insertions(+), 641 deletions(-) diff --git a/NEWS b/NEWS index 63f4a50..992a1be 100644 --- a/NEWS +++ b/NEWS @@ -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 => 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). diff --git a/include/GArf.hpp b/include/GArf.hpp index 9d30270..af450a6 100644 --- a/include/GArf.hpp +++ b/include/GArf.hpp @@ -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; diff --git a/include/GEnergies.hpp b/include/GEnergies.hpp index 9141430..d06d318 100644 --- a/include/GEnergies.hpp +++ b/include/GEnergies.hpp @@ -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; diff --git a/include/GHealpix.hpp b/include/GHealpix.hpp index ec2504c..ba6e961 100644 --- a/include/GHealpix.hpp +++ b/include/GHealpix.hpp @@ -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; diff --git a/include/GPha.hpp b/include/GPha.hpp index 82a094d..28ce207 100644 --- a/include/GPha.hpp +++ b/include/GPha.hpp @@ -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; diff --git a/include/GPhotons.hpp b/include/GPhotons.hpp index a145c83..61ae4bc 100644 --- a/include/GPhotons.hpp +++ b/include/GPhotons.hpp @@ -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 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 */ diff --git a/include/GRmf.hpp b/include/GRmf.hpp index 6f9b060..d6aef73 100644 --- a/include/GRmf.hpp +++ b/include/GRmf.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; diff --git a/include/GSkyProjection.hpp b/include/GSkyProjection.hpp index aac4f56..26875d8 100644 --- a/include/GSkyProjection.hpp +++ b/include/GSkyProjection.hpp @@ -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; diff --git a/include/GSkymap.hpp b/include/GSkymap.hpp index fbfbe78..5f85cd1 100644 --- a/include/GSkymap.hpp +++ b/include/GSkymap.hpp @@ -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; diff --git a/include/GTimes.hpp b/include/GTimes.hpp index 0c35893..3e50db4 100644 --- a/include/GTimes.hpp +++ b/include/GTimes.hpp @@ -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 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 */ diff --git a/include/GWcs.hpp b/include/GWcs.hpp index 399c151..4d2349d 100644 --- a/include/GWcs.hpp +++ b/include/GWcs.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; diff --git a/inst/com/src/GCOMEventCube.cpp b/inst/com/src/GCOMEventCube.cpp index 01bca0f..066a8da 100644 --- a/inst/com/src/GCOMEventCube.cpp +++ b/inst/com/src/GCOMEventCube.cpp @@ -417,11 +417,11 @@ void GCOMEventCube::read(const GFits& 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& 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& 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& 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& fits) void GCOMEventCube::write(GFits& file) const { // Write cube - m_map.write(&file); + m_map.write(file); // Write energy boundaries diff --git a/inst/com/src/GCOMObservation.cpp b/inst/com/src/GCOMObservation.cpp index 1646a8d..8f1b3d2 100644 --- a/inst/com/src/GCOMObservation.cpp +++ b/inst/com/src/GCOMObservation.cpp @@ -734,11 +734,11 @@ void GCOMObservation::load_drb(const std::string& 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& 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& 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); diff --git a/inst/cta/pyext/GCTAInstDir.i b/inst/cta/pyext/GCTAInstDir.i index 7f352f9..6201cab 100644 --- a/inst/cta/pyext/GCTAInstDir.i +++ b/inst/cta/pyext/GCTAInstDir.i @@ -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); }; diff --git a/inst/cta/src/GCTAEventCube.cpp b/inst/cta/src/GCTAEventCube.cpp index 8d40ab6..04d2c51 100644 --- a/inst/cta/src/GCTAEventCube.cpp +++ b/inst/cta/src/GCTAEventCube.cpp @@ -396,7 +396,7 @@ void GCTAEventCube::read(const GFits& 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(); diff --git a/inst/lat/src/GLATEventCube.cpp b/inst/lat/src/GLATEventCube.cpp index c7e634a..e0ecfb7 100644 --- a/inst/lat/src/GLATEventCube.cpp +++ b/inst/lat/src/GLATEventCube.cpp @@ -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 diff --git a/inst/lat/src/GLATLtCubeMap.cpp b/inst/lat/src/GLATLtCubeMap.cpp index 702649e..33f36e6 100644 --- a/inst/lat/src/GLATLtCubeMap.cpp +++ b/inst/lat/src/GLATLtCubeMap.cpp @@ -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"))); diff --git a/pyext/GArf.i b/pyext/GArf.i index d3ce275..fc9377c 100644 --- a/pyext/GArf.i +++ b/pyext/GArf.i @@ -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; }; diff --git a/pyext/GEnergies.i b/pyext/GEnergies.i index b694315..e056fa6 100644 --- a/pyext/GEnergies.i +++ b/pyext/GEnergies.i @@ -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; }; diff --git a/pyext/GHealpix.i b/pyext/GHealpix.i index 5c231a0..e39ebe3 100644 --- a/pyext/GHealpix.i +++ b/pyext/GHealpix.i @@ -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; diff --git a/pyext/GPha.i b/pyext/GPha.i index e29f396..dcb910a 100644 --- a/pyext/GPha.i +++ b/pyext/GPha.i @@ -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; }; diff --git a/pyext/GRmf.i b/pyext/GRmf.i index 9524e88..ccb045f 100644 --- a/pyext/GRmf.i +++ b/pyext/GRmf.i @@ -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; }; diff --git a/pyext/GSkyProjection.i b/pyext/GSkyProjection.i index 40876c8..46c006c 100644 --- a/pyext/GSkyProjection.i +++ b/pyext/GSkyProjection.i @@ -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; diff --git a/pyext/GSkymap.i b/pyext/GSkymap.i index 6892900..5cc6f90 100644 --- a/pyext/GSkymap.i +++ b/pyext/GSkymap.i @@ -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; }; diff --git a/pyext/GWcs.i b/pyext/GWcs.i index 463738e..2f63c29 100644 --- a/pyext/GWcs.i +++ b/pyext/GWcs.i @@ -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; diff --git a/src/obs/GEnergies.cpp b/src/obs/GEnergies.cpp index 1b2602d..e67a475 100644 --- a/src/obs/GEnergies.cpp +++ b/src/obs/GEnergies.cpp @@ -363,11 +363,11 @@ void GEnergies::load(const std::string& filename, const std::string& 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& 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& 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& extname) const table->extname(extname); // Write to FITS file - file->append(*table); + file.append(*table); // Free table delete table; diff --git a/src/sky/GHealpix.cpp b/src/sky/GHealpix.cpp index 2c54313..066c5de 100644 --- a/src/sky/GHealpix.cpp +++ b/src/sky/GHealpix.cpp @@ -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& 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; diff --git a/src/sky/GSkymap.cpp b/src/sky/GSkymap.cpp index d6e8c5c..071bf78 100644 --- a/src/sky/GSkymap.cpp +++ b/src/sky/GSkymap.cpp @@ -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& 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(hdu)); + if (hdu.hascard("PIXTYPE") && hdu.string("PIXTYPE") == "HEALPIX") { + read_healpix(static_cast(hdu)); loaded = true; break; } @@ -963,17 +962,17 @@ void GSkymap::load(const std::string& 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(hdu)); + read_wcs(static_cast(hdu)); loaded = true; break; @@ -1035,39 +1034,32 @@ void GSkymap::save(const std::string& 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(hdu)); + loaded = true; + } - // If PIXTYPE keyword equals "HEALPIX" then load map - if (hdu->hascard("PIXTYPE") && hdu->string("PIXTYPE") == "HEALPIX") { - read_healpix(static_cast(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(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(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& wcs, const std::string& 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& wcs, const std::string& 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(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(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 diff --git a/src/sky/GWcs.cpp b/src/sky/GWcs.cpp index 4378c63..4e8e1d7 100644 --- a/src/sky/GWcs.cpp +++ b/src/sky/GWcs.cpp @@ -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& 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& GWcs::operator=(const GWcs& 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; diff --git a/src/xspec/GArf.cpp b/src/xspec/GArf.cpp index 303ec5c..e371920 100644 --- a/src/xspec/GArf.cpp +++ b/src/xspec/GArf.cpp @@ -259,7 +259,7 @@ void GArf::load(const std::string& 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& filename, const bool& 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; diff --git a/src/xspec/GPha.cpp b/src/xspec/GPha.cpp index e1b45b4..761991b 100644 --- a/src/xspec/GPha.cpp +++ b/src/xspec/GPha.cpp @@ -339,7 +339,7 @@ void GPha::load(const std::string& 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& filename, const bool& 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; diff --git a/src/xspec/GRmf.cpp b/src/xspec/GRmf.cpp index e4a29c7..84ff777 100644 --- a/src/xspec/GRmf.cpp +++ b/src/xspec/GRmf.cpp @@ -276,7 +276,7 @@ void GRmf::load(const std::string& 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& filename, const bool& 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; -- libgit2 0.21.4