diff --git a/configure.ac b/configure.ac index 92eb899..0332901 100644 --- a/configure.ac +++ b/configure.ac @@ -659,8 +659,8 @@ if test "x$openmp_ok" = "xyes" ; then # dynamically if test "$OPENMP_CXXFLAGS" = "-fopenmp" -a "$ac_sys_system" = "Linux"; then AC_MSG_CHECKING(whether we can dynamically load libgomp) - SAVED_LDFLAGS=$LDFLAGS - LDFLAGS="$LDFLAGS -ldl" + SAVED_LIBS=$LIBS + LIBS="-ldl" AC_TRY_RUN(AC_LANG_PROGRAM([#include ], [if (dlopen("libgomp.so.1", RTLD_NOW) == 0) return 1;]), [AC_MSG_RESULT(yes)], @@ -668,7 +668,7 @@ if test "x$openmp_ok" = "xyes" ; then AC_MSG_RESULT(no) AC_MSG_WARN([OpenMP detected but shared libgomp.so.1 library cannot be dlopen()ed.]) AC_MSG_WARN([If you use a gcc 4.2.x compiler, this is a known feature. Try to use gcc >= 4.3.])]) - LDFLAGS=$SAVED_LDFLAGS + LIBS=$SAVED_LIBS fi # Special Darwin test to see whether OpenMP compiles without adding the diff --git a/include/GEbounds.hpp b/include/GEbounds.hpp index 67300e4..e2917c2 100644 --- a/include/GEbounds.hpp +++ b/include/GEbounds.hpp @@ -76,7 +76,7 @@ public: GEnergy emax(int inx) const; GEnergy emean(int inx) const; GEnergy elogmean(int inx) const; - bool isin(const GEnergy& eng) const; + bool contains(const GEnergy& eng) const; std::string print(void) const; diff --git a/inst/cta/include/GCTAEventBin.hpp b/inst/cta/include/GCTAEventBin.hpp index c59d240..29d145f 100644 --- a/inst/cta/include/GCTAEventBin.hpp +++ b/inst/cta/include/GCTAEventBin.hpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GCTAEventBin.hpp - CTA event bin class * + * GCTAEventBin.hpp - CTA event bin class * * ----------------------------------------------------------------------- * - * copyright (C) 2010 by Jurgen Knodlseder * + * copyright (C) 2010-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,8 +20,8 @@ ***************************************************************************/ /** * @file GCTAEventBin.hpp - * @brief GCTAEventBin class interface definition. - * @author J. Knodlseder + * @brief CTA event bin class interface definition + * @author Juergen Knoedlseder */ #ifndef GCTAEVENTBIN_HPP @@ -67,18 +67,18 @@ public: virtual void clear(void); virtual GCTAEventBin* clone(void) const; virtual double size(void) const; - virtual const GCTAInstDir& dir(void) const { return *m_dir; } - virtual const GEnergy& energy(void) const { return *m_energy; } - virtual const GTime& time(void) const { return *m_time; } - virtual double counts(void) const { return *m_counts; } + virtual const GCTAInstDir& dir(void) const; + virtual const GEnergy& energy(void) const; + virtual const GTime& time(void) const; + virtual double counts(void) const; virtual double error(void) const; - virtual void counts(const double& counts) { *m_counts=counts; } + virtual void counts(const double& counts); virtual std::string print(void) const; // Other methods - const double& omega(void) const { return *m_omega; } - const GEnergy& ewidth(void) const { return *m_ewidth; } - const double& ontime(void) const { return *m_ontime; } + const double& omega(void) const; + const GEnergy& ewidth(void) const; + const double& ontime(void) const; protected: // Protected methods diff --git a/inst/cta/include/GCTAException.hpp b/inst/cta/include/GCTAException.hpp index 0a0047f..b5ed5a9 100644 --- a/inst/cta/include/GCTAException.hpp +++ b/inst/cta/include/GCTAException.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAException.hpp - CTA exception handler * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2011 by Jurgen Knodlseder * + * copyright (C) 2010-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file GCTAException.hpp * @brief CTA exception handler interface definition. - * @author J. Knodlseder + * @author Juergen Knoedlseder */ #ifndef GCTAEXCEPTION_HPP @@ -65,6 +65,12 @@ public: std::string message = ""); }; + // Event bin exceptions + class no_member : public GCTAExceptionHandler { + public: + no_member(const std::string& origin, const std::string& message); + }; + // Observation exceptions class no_response : public GCTAExceptionHandler { public: diff --git a/inst/cta/src/GCTAEventBin.cpp b/inst/cta/src/GCTAEventBin.cpp index 7f8be1a..e12267d 100644 --- a/inst/cta/src/GCTAEventBin.cpp +++ b/inst/cta/src/GCTAEventBin.cpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GCTAEventBin.cpp - CTA event bin class * + * GCTAEventBin.cpp - CTA event bin class * * ----------------------------------------------------------------------- * - * copyright (C) 2010 by Jurgen Knodlseder * + * copyright (C) 2010-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,8 +20,8 @@ ***************************************************************************/ /** * @file GCTAEventBin.cpp - * @brief GCTAEventBin class implementation. - * @author J. Knodlseder + * @brief CTA event bin class implementation + * @author Juergen Knoedlseder */ /* __ Includes ___________________________________________________________ */ @@ -31,9 +31,18 @@ #include #include #include "GCTAEventBin.hpp" +#include "GCTAException.hpp" #include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ +#define G_DIR "GCTAEventBin::dir()" +#define G_ENERGY "GCTAEventBin::energy()" +#define G_TIME "GCTAEventBin::time()" +#define G_COUNTS_GET "GCTAEventBin::counts()" +#define G_COUNTS_SET "GCTAEventBin::counts(double&)" +#define G_OMEGA "GCTAEventBin::omega()" +#define G_EWIDTH "GCTAEventBin::ewidth()" +#define G_ONTIME "GCTAEventBin::ontime()" /* __ Macros _____________________________________________________________ */ @@ -185,8 +194,128 @@ double GCTAEventBin::size(void) const /***********************************************************************//** + * @brief Return instrument direction of event bin + * + * @return Instrument direction of event bin + * + * @exception GCTAException::no_member + * Invalid instrument direction pointer. + * + * Returns reference to the instrument direction of the event bin. + ***************************************************************************/ +const GCTAInstDir& GCTAEventBin::dir(void) const +{ + // Throw an exception if instrument direction pointer is not valid + if (m_dir == NULL) { + throw GCTAException::no_member(G_DIR, + "Invalid instrument direction pointer."); + } + + // Return instrument direction + return *m_dir; +} + + +/***********************************************************************//** + * @brief Return energy of event bin + * + * @return Energy of event bin + * + * @exception GCTAException::no_member + * Invalid energy pointer. + * + * Returns reference to the energy of the event bin. + ***************************************************************************/ +const GEnergy& GCTAEventBin::energy(void) const +{ + // Throw an exception if energy pointer is not valid + if (m_energy == NULL) { + throw GCTAException::no_member(G_ENERGY, + "Invalid energy pointer."); + } + + // Return energy + return *m_energy; +} + + +/***********************************************************************//** + * @brief Return time of event bin + * + * @return Time of event bin + * + * @exception GCTAException::no_member + * Invalid time pointer. + * + * Returns reference to the time of the event bin. + ***************************************************************************/ +const GTime& GCTAEventBin::time(void) const +{ + // Throw an exception if time pointer is not valid + if (m_time == NULL) { + throw GCTAException::no_member(G_TIME, + "Invalid time pointer."); + } + + // Return time + return *m_time; +} + + +/***********************************************************************//** + * @brief Return number of counts in event bin + * + * @return Number of counts in event bin + * + * @exception GCTAException::no_member + * Invalid counts pointer. + * + * Returns reference to the number of counts in the event bin. + ***************************************************************************/ +double GCTAEventBin::counts(void) const +{ + // Throw an exception if counts pointer is not valid + if (m_counts == NULL) { + throw GCTAException::no_member(G_COUNTS_GET, + "Invalid counts pointer."); + } + + // Return counts + return *m_counts; +} + + +/***********************************************************************//** + * @brief Set number of counts in event bin + * + * @param[in] counts Number of counts. + * + * @exception GCTAException::no_member + * Invalid counts pointer. + * + * Set the number of counts in the event bin. + ***************************************************************************/ +void GCTAEventBin::counts(const double& counts) +{ + // Throw an exception if counts pointer is not valid + if (m_counts == NULL) { + throw GCTAException::no_member(G_COUNTS_SET, + "Invalid counts pointer."); + } + + // Set number of counts in event bin + *m_counts = counts; + + // Return + return; +} + + +/***********************************************************************//** * @brief Return error in number of counts * + * @return Error in number of counts in event bin + * * Returns \f$\sqrt(counts+delta)\f$ as the uncertainty in the number of * counts in the bin. Adding delta avoids uncertainties of 0 which will * lead in the optimisation step to the exlusion of the corresponding bin. @@ -206,6 +335,75 @@ double GCTAEventBin::error(void) const /***********************************************************************//** + * @brief Return solid angle of event bin + * + * @return Solid angle of event bin + * + * @exception GCTAException::no_member + * Invalid solid angle pointer. + * + * Returns reference to the solid angle of the event bin. + ***************************************************************************/ +const double& GCTAEventBin::omega(void) const +{ + // Throw an exception if solid angle pointer is not valid + if (m_omega == NULL) { + throw GCTAException::no_member(G_OMEGA, + "Invalid solid angle pointer."); + } + + // Return solid angle + return *m_omega; +} + + +/***********************************************************************//** + * @brief Return energy width of event bin + * + * @return Energy width of event bin + * + * @exception GCTAException::no_member + * Invalid energy width pointer. + * + * Returns reference to the energy width of the event bin. + ***************************************************************************/ +const GEnergy& GCTAEventBin::ewidth(void) const +{ + // Throw an exception if energy width pointer is not valid + if (m_ewidth == NULL) { + throw GCTAException::no_member(G_EWIDTH, + "Invalid energy width pointer."); + } + + // Return energy width + return *m_ewidth; +} + + +/***********************************************************************//** + * @brief Return ontime of event bin + * + * @return Ontime of event bin + * + * @exception GCTAException::no_member + * Invalid ontime pointer. + * + * Returns reference to the ontime of the event bin. + ***************************************************************************/ +const double& GCTAEventBin::ontime(void) const +{ + // Throw an exception if ontime pointer is not valid + if (m_ontime == NULL) { + throw GCTAException::no_member(G_ONTIME, + "Invalid ontime pointer."); + } + + // Return ontime + return *m_ontime; +} + + +/***********************************************************************//** * @brief Print event information ***************************************************************************/ std::string GCTAEventBin::print(void) const diff --git a/inst/cta/src/GCTAException.cpp b/inst/cta/src/GCTAException.cpp index 897c696..feca454 100644 --- a/inst/cta/src/GCTAException.cpp +++ b/inst/cta/src/GCTAException.cpp @@ -61,6 +61,26 @@ GCTAException::file_open_error::file_open_error(std::string origin, /***********************************************************************//** + * @brief Member not set. + * + * @param[in] origin Name of method that has thrown the exception. + * @param[in] message Error message. + ***************************************************************************/ +GCTAException::no_member::no_member(const std::string& origin, + const std::string& message) +{ + // Set origin and message + m_origin = origin; + m_message = message; + + // Return + return; +} + + + + +/***********************************************************************//** * @brief Pointing not set. * * @param[in] origin Name of method that has thrown the exception. diff --git a/inst/lat/include/GLATEventBin.hpp b/inst/lat/include/GLATEventBin.hpp index 6d9aaa9..15459a3 100644 --- a/inst/lat/include/GLATEventBin.hpp +++ b/inst/lat/include/GLATEventBin.hpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GLATEventBin.hpp - LAT event bin class * + * GLATEventBin.hpp - LAT event bin class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2011 by Jurgen Knodlseder * + * copyright (C) 2009-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file GLATEventBin.hpp * @brief Fermi-LAT event bin class interface definition - * @author J. Knodlseder + * @author Juergen Knoedlseder */ #ifndef GLATEVENTBIN_HPP @@ -62,18 +62,18 @@ public: virtual void clear(void); virtual GLATEventBin* clone(void) const; virtual double size(void) const; - virtual const GLATInstDir& dir(void) const { return *m_dir; } - virtual const GEnergy& energy(void) const { return *m_energy; } - virtual const GTime& time(void) const { return *m_time; } - virtual double counts(void) const { return *m_counts; } + virtual const GLATInstDir& dir(void) const; + virtual const GEnergy& energy(void) const; + virtual const GTime& time(void) const; + virtual double counts(void) const; virtual double error(void) const; - virtual void counts(const double& counts) { *m_counts=counts; } + virtual void counts(const double& counts); virtual std::string print(void) const; // Other methods - const double& omega(void) const { return *m_omega; } - const GEnergy& ewidth(void) const { return *m_ewidth; } - const double& ontime(void) const { return *m_ontime; } + const double& omega(void) const; + const GEnergy& ewidth(void) const; + const double& ontime(void) const; const int& index(void) const { return m_index; } const int& ipix(void) const { return m_ipix; } const int& ieng(void) const { return m_ieng; } diff --git a/inst/lat/include/GLATException.hpp b/inst/lat/include/GLATException.hpp index 798e0ef..87141af 100644 --- a/inst/lat/include/GLATException.hpp +++ b/inst/lat/include/GLATException.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GLATException.hpp - LAT exception handler * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2011 by Jurgen Knodlseder * + * copyright (C) 2010-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file GLATException.hpp * @brief LAT exception handler interface definition. - * @author J. Knodlseder + * @author Juergen Knoedlseder */ #ifndef GLATEXCEPTION_HPP @@ -64,6 +64,10 @@ public: file_open_error(std::string origin, std::string filename, std::string message = ""); }; + class no_member : public GLATExceptionHandler { + public: + no_member(const std::string& origin, const std::string& message); + }; class no_sky : public GLATExceptionHandler { public: no_sky(std::string origin, std::string message = ""); diff --git a/inst/lat/src/GLATAeff.cpp b/inst/lat/src/GLATAeff.cpp index 69b3303..247d0dc 100644 --- a/inst/lat/src/GLATAeff.cpp +++ b/inst/lat/src/GLATAeff.cpp @@ -539,8 +539,8 @@ void GLATAeff::copy_members(const GLATAeff& aeff) m_back = aeff.m_back; // Clone functors - m_eff_func1 = aeff.m_eff_func1->clone(); - m_eff_func2 = aeff.m_eff_func2->clone(); + m_eff_func1 = (aeff.m_eff_func1 != NULL) ? aeff.m_eff_func1->clone() : NULL; + m_eff_func2 = (aeff.m_eff_func2 != NULL) ? aeff.m_eff_func2->clone() : NULL; // Return return; diff --git a/inst/lat/src/GLATEventBin.cpp b/inst/lat/src/GLATEventBin.cpp index 3bbadd7..f5dbca6 100644 --- a/inst/lat/src/GLATEventBin.cpp +++ b/inst/lat/src/GLATEventBin.cpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GLATEventBin.cpp - LAT event bin class * + * GLATEventBin.cpp - LAT event bin class * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2010 by Jurgen Knodlseder * + * copyright (C) 2009-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -20,8 +20,8 @@ ***************************************************************************/ /** * @file GLATEventBin.cpp - * @brief GLATEventBin class implementation. - * @author J. Knodlseder + * @brief Fermi-LAT event bin class implementation + * @author Juergen Knodlseder */ /* __ Includes ___________________________________________________________ */ @@ -31,9 +31,18 @@ #include #include #include "GLATEventBin.hpp" +#include "GLATException.hpp" #include "GTools.hpp" /* __ Method name definitions ____________________________________________ */ +#define G_DIR "GLATEventBin::dir()" +#define G_ENERGY "GLATEventBin::energy()" +#define G_TIME "GLATEventBin::time()" +#define G_COUNTS_GET "GLATEventBin::counts()" +#define G_COUNTS_SET "GLATEventBin::counts(double&)" +#define G_OMEGA "GLATEventBin::omega()" +#define G_EWIDTH "GLATEventBin::ewidth()" +#define G_ONTIME "GLATEventBin::ontime()" /* __ Macros _____________________________________________________________ */ @@ -101,7 +110,8 @@ GLATEventBin::~GLATEventBin(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] bin Event bin which should be assigned. + * @param[in] bin LAT event bin. + * @return Lat event bin. ***************************************************************************/ GLATEventBin& GLATEventBin::operator= (const GLATEventBin& bin) { @@ -157,7 +167,9 @@ void GLATEventBin::clear(void) /***********************************************************************//** * @brief Clone instance -***************************************************************************/ + * + * @return Deep copy of LAT event bin. + ***************************************************************************/ GLATEventBin* GLATEventBin::clone(void) const { return new GLATEventBin(*this); @@ -167,6 +179,8 @@ GLATEventBin* GLATEventBin::clone(void) const /***********************************************************************//** * @brief Return size of event bin * + * @return Size of event bin in units of sr MeV s. + * * The size of the event bin (units sr MeV s) is given by * \f[size = \Omega \times \Delta E \times \Delta T\f] * where @@ -185,8 +199,128 @@ double GLATEventBin::size(void) const /***********************************************************************//** + * @brief Return instrument direction of event bin + * + * @return Instrument direction of event bin. + * + * @exception GLATException::no_member + * Invalid instrument direction pointer. + * + * Returns reference to the instrument direction of the event bin. + ***************************************************************************/ +const GLATInstDir& GLATEventBin::dir(void) const +{ + // Throw an exception if instrument direction pointer is not valid + if (m_dir == NULL) { + throw GLATException::no_member(G_DIR, + "Invalid instrument direction pointer."); + } + + // Return instrument direction + return *m_dir; +} + + +/***********************************************************************//** + * @brief Return energy of event bin + * + * @return Energy of event bin. + * + * @exception GLATException::no_member + * Invalid energy pointer. + * + * Returns reference to the energy of the event bin. + ***************************************************************************/ +const GEnergy& GLATEventBin::energy(void) const +{ + // Throw an exception if energy pointer is not valid + if (m_energy == NULL) { + throw GLATException::no_member(G_ENERGY, + "Invalid energy pointer."); + } + + // Return energy + return *m_energy; +} + + +/***********************************************************************//** + * @brief Return time of event bin + * + * @return Time of event bin. + * + * @exception GLATException::no_member + * Invalid time pointer. + * + * Returns reference to the time of the event bin. + ***************************************************************************/ +const GTime& GLATEventBin::time(void) const +{ + // Throw an exception if time pointer is not valid + if (m_time == NULL) { + throw GLATException::no_member(G_TIME, + "Invalid time pointer."); + } + + // Return time + return *m_time; +} + + +/***********************************************************************//** + * @brief Return number of counts in event bin + * + * @return Number of counts in event bin. + * + * @exception GLATException::no_member + * Invalid counts pointer. + * + * Returns reference to the number of counts in the event bin. + ***************************************************************************/ +double GLATEventBin::counts(void) const +{ + // Throw an exception if counts pointer is not valid + if (m_counts == NULL) { + throw GLATException::no_member(G_COUNTS_GET, + "Invalid counts pointer."); + } + + // Return counts + return *m_counts; +} + + +/***********************************************************************//** + * @brief Set number of counts in event bin + * + * @param[in] counts Number of counts. + * + * @exception GLATException::no_member + * Invalid counts pointer. + * + * Set the number of counts in the event bin. + ***************************************************************************/ +void GLATEventBin::counts(const double& counts) +{ + // Throw an exception if counts pointer is not valid + if (m_counts == NULL) { + throw GLATException::no_member(G_COUNTS_SET, + "Invalid counts pointer."); + } + + // Set number of counts in event bin + *m_counts = counts; + + // Return + return; +} + + +/***********************************************************************//** * @brief Return error in number of counts * + * @return Error in number of counts in event bin. + * * Returns \f$\sqrt(counts+delta)\f$ as the uncertainty in the number of * counts in the bin. Adding delta avoids uncertainties of 0 which will * lead in the optimisation step to the exlusion of the corresponding bin. @@ -206,7 +340,78 @@ double GLATEventBin::error(void) const /***********************************************************************//** + * @brief Return solid angle of event bin + * + * @return Solid angle of event bin. + * + * @exception GLATException::no_member + * Invalid solid angle pointer. + * + * Returns reference to the solid angle of the event bin. + ***************************************************************************/ +const double& GLATEventBin::omega(void) const +{ + // Throw an exception if solid angle pointer is not valid + if (m_omega == NULL) { + throw GLATException::no_member(G_OMEGA, + "Invalid solid angle pointer."); + } + + // Return solid angle + return *m_omega; +} + + +/***********************************************************************//** + * @brief Return energy width of event bin + * + * @return Energy width of event bin. + * + * @exception GLATException::no_member + * Invalid energy width pointer. + * + * Returns reference to the energy width of the event bin. + ***************************************************************************/ +const GEnergy& GLATEventBin::ewidth(void) const +{ + // Throw an exception if energy width pointer is not valid + if (m_ewidth == NULL) { + throw GLATException::no_member(G_EWIDTH, + "Invalid energy width pointer."); + } + + // Return energy width + return *m_ewidth; +} + + +/***********************************************************************//** + * @brief Return ontime of event bin + * + * @return Ontime of event bin. + * + * @exception GLATException::no_member + * Invalid ontime pointer. + * + * Returns reference to the ontime of the event bin. + ***************************************************************************/ +const double& GLATEventBin::ontime(void) const +{ + // Throw an exception if ontime pointer is not valid + if (m_ontime == NULL) { + throw GLATException::no_member(G_ONTIME, + "Invalid ontime pointer."); + } + + // Return ontime + return *m_ontime; +} + + +/***********************************************************************//** * @brief Print event information + * + * @return String containing number of counts in event bin. ***************************************************************************/ std::string GLATEventBin::print(void) const { @@ -283,10 +488,3 @@ void GLATEventBin::free_members(void) // Return return; } - - -/*========================================================================== - = = - = Friends = - = = - ==========================================================================*/ diff --git a/inst/lat/src/GLATEventCube.cpp b/inst/lat/src/GLATEventCube.cpp index 1dde788..a009447 100644 --- a/inst/lat/src/GLATEventCube.cpp +++ b/inst/lat/src/GLATEventCube.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * GLATEventCube.cpp - LAT event cube class * + * GLATEventCube.cpp - LAT event cube class * * ----------------------------------------------------------------------- * * copyright (C) 2009-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -21,7 +21,7 @@ /** * @file GLATEventCube.cpp * @brief LAT event cube class implementation - * @author J. Knoedlseder + * @author Juergen Knoedlseder */ /* __ Includes ___________________________________________________________ */ @@ -111,6 +111,7 @@ GLATEventCube::~GLATEventCube(void) * @brief Assignment operator * * @param[in] cube LAT event cube. + * @return LAT event cube. ***************************************************************************/ GLATEventCube& GLATEventCube::operator= (const GLATEventCube& cube) { @@ -140,6 +141,7 @@ GLATEventCube& GLATEventCube::operator= (const GLATEventCube& cube) * @brief Event bin access operator * * @param[in] index Event index [0,...,size()-1]. + * @return Pointer to event bin. * * Returns pointer to an event bin. ***************************************************************************/ @@ -157,6 +159,7 @@ GLATEventBin* GLATEventCube::operator[](const int& index) * @brief Event bin access operator (const version) * * @param[in] index Event index [0,...,size()-1]. + * @return Pointer to event bin. * * Returns pointer to an event bin. ***************************************************************************/ @@ -200,6 +203,8 @@ void GLATEventCube::clear(void) /***********************************************************************//** * @brief Clone instance + * + * @return Deep copy of LAT event cube. ***************************************************************************/ GLATEventCube* GLATEventCube::clone(void) const { @@ -209,6 +214,8 @@ GLATEventCube* GLATEventCube::clone(void) const /***********************************************************************//** * @brief Return number of bins in event cube + * + * @return Number of bins in event cube. ***************************************************************************/ int GLATEventCube::size(void) const { @@ -222,6 +229,8 @@ int GLATEventCube::size(void) const /***********************************************************************//** * @brief Return dimension of event cube + * + * @return Number of dimensions in event cube (either 2 or 3). ***************************************************************************/ int GLATEventCube::dim(void) const { @@ -236,7 +245,8 @@ int GLATEventCube::dim(void) const /***********************************************************************//** * @brief Return number of bins in axis * - * @param[in] axis Axis. + * @param[in] axis Axis [0,...,dim()-1] + * @return Number of bins in specified axis * * @exception GException::out_of_range * Axis is out of range. @@ -247,8 +257,9 @@ int GLATEventCube::naxis(int axis) const { // Optionally check if the axis is valid #if defined(G_RANGE_CHECK) - if (axis < 0 || axis >= dim()) + if (axis < 0 || axis >= dim()) { throw GException::out_of_range(G_NAXIS, axis, 0, dim()-1); + } #endif // Set result @@ -368,6 +379,8 @@ void GLATEventCube::write(GFits& file) const /***********************************************************************//** * @brief Return number of events in cube + * + * @return Total number of events in event cube. ***************************************************************************/ int GLATEventCube::number(void) const { @@ -379,8 +392,9 @@ int GLATEventCube::number(void) const // Sum event cube if (size() > 0 && pixels != NULL) { - for (int i = 0; i < size(); ++i) + for (int i = 0; i < size(); ++i) { number += pixels[i]; + } } // Return @@ -406,6 +420,8 @@ void GLATEventCube::map(const GSkymap& map) /***********************************************************************//** * @brief Print event cube information + * + * @return String containing event cube information. ***************************************************************************/ std::string GLATEventCube::print(void) const { @@ -419,17 +435,29 @@ std::string GLATEventCube::print(void) const result.append(str(m_map.nx())+" x "+str(m_map.ny())); result.append("\n"+parformat("Number of energy bins")+str(ebins())); result.append("\n"+parformat("Number of events")+str(number())); + + // Append time interval result.append("\n"+parformat("Time interval")); - if (gti().size() > 0) + if (gti().size() > 0) { result.append(str(tstart().met())+" - "+str(tstop().met())); - else + } + else { result.append("not defined"); + } + + // Append energy range result.append("\n"+parformat("Energy range")); - if (ebounds().size() > 0) + if (ebounds().size() > 0) { result.append(emin().print()+" - "+emax().print()); - else + } + else { result.append("not defined"); - result.append("\n"+m_map.wcs()->print()); + } + + // Append WCS + if (m_map.wcs() != NULL) { + result.append("\n"+m_map.wcs()->print()); + } // Append source maps result.append("\n"+parformat("Number of source maps")+str(m_srcmap.size())); @@ -450,7 +478,8 @@ std::string GLATEventCube::print(void) const /***********************************************************************//** * @brief Return name of diffuse model * - * @param[in] index Diffuse model index (starting from 0). + * @param[in] index Diffuse model index [0,...,ndiffrsp()-1]. + * @return Name of diffuse model. * * @exception GException::out_of_range * Model index out of valid range. @@ -461,8 +490,9 @@ std::string GLATEventCube::diffname(const int& index) const { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (index < 0 || index >= ndiffrsp()) + if (index < 0 || index >= ndiffrsp()) { throw GException::out_of_range(G_DIFFNAME, index, 0, ndiffrsp()-1); + } #endif // Return @@ -473,7 +503,8 @@ std::string GLATEventCube::diffname(const int& index) const /***********************************************************************//** * @brief Return diffuse response map * - * @param[in] index Diffuse model index (starting from 0). + * @param[in] index Diffuse model index [0,...,ndiffrsp()-1]. + * @return Pointer to diffuse response map. * * @exception GException::out_of_range * Model index out of valid range. @@ -484,8 +515,9 @@ GSkymap* GLATEventCube::diffrsp(const int& index) const { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (index < 0 || index >= ndiffrsp()) + if (index < 0 || index >= ndiffrsp()) { throw GException::out_of_range(G_DIFFRSP, index, 0, ndiffrsp()-1); + } #endif // Return @@ -498,6 +530,7 @@ GSkymap* GLATEventCube::diffrsp(const int& index) const * direction that fits spatially into the event cube * * @param[in] srcDir Source direction. + * @return Maximum radius in degrees that fully fits into event cube. * * By computing the sky directions of the event cube boundaries, the maximum * radius is computed that fits fully within the event cube. This method is @@ -513,8 +546,9 @@ double GLATEventCube::maxrad(const GSkyDir& srcDir) const for (int ix = 0; ix < nx(); ++ix) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); double distance = m_map.xy2dir(pixel).dist_deg(srcDir); - if (distance < radius) + if (distance < radius) { radius = distance; + } } // Move along lower edge in longitude @@ -522,8 +556,9 @@ double GLATEventCube::maxrad(const GSkyDir& srcDir) const for (int ix = 0; ix < nx(); ++ix) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); double distance = m_map.xy2dir(pixel).dist_deg(srcDir); - if (distance < radius) + if (distance < radius) { radius = distance; + } } // Move along left edge in latitude @@ -531,8 +566,9 @@ double GLATEventCube::maxrad(const GSkyDir& srcDir) const for (int iy = 0; iy < ny(); ++iy) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); double distance = m_map.xy2dir(pixel).dist_deg(srcDir); - if (distance < radius) + if (distance < radius) { radius = distance; + } } // Move along right edge in latitude @@ -540,8 +576,9 @@ double GLATEventCube::maxrad(const GSkyDir& srcDir) const for (int iy = 0; iy < ny(); ++iy) { GSkyPixel pixel = GSkyPixel(double(ix), double(iy)); double distance = m_map.xy2dir(pixel).dist_deg(srcDir); - if (distance < radius) + if (distance < radius) { radius = distance; + } } // Return radius @@ -664,18 +701,21 @@ void GLATEventCube::read_srcmap(const GFitsImage* hdu) map->read(hdu); // Check that source map WCS is consistent with counts map WCS - if (*(m_map.wcs()) != *(map->wcs())) + if (*(m_map.wcs()) != *(map->wcs())) { throw GLATException::wcs_incompatible(G_READ_SRCMAP, hdu->extname()); + } // Check that source map dimension is consistent with counts map // dimension if (m_map.nx() != map->nx() || - m_map.ny() != map->ny()) + m_map.ny() != map->ny()) { throw GLATException::wcs_incompatible(G_READ_SRCMAP, hdu->extname()); + } // Check that source map has required number of energy bins - if (m_map.nmaps()+1 != map->nmaps()) + if (m_map.nmaps()+1 != map->nmaps()) { throw GLATException::wcs_incompatible(G_READ_SRCMAP, hdu->extname()); + } // Append source map to list of maps m_srcmap.push_back(map); @@ -837,9 +877,10 @@ void GLATEventCube::set_energies(void) void GLATEventCube::set_times(void) { // Throw an error if GTI is empty - if (m_gti.size() < 1) + if (m_gti.size() < 1) { throw GLATException::no_gti(G_SET_TIMES, "Every LAT event cube needs" " associated GTIs to allow the computation of the ontime."); + } // Compute mean time m_time = 0.5 * (gti().tstart() + gti().tstop()); @@ -874,17 +915,20 @@ void GLATEventCube::set_bin(const int& index) { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (index < 0 || index >= size()) + if (index < 0 || index >= size()) { throw GException::out_of_range(G_SET_BIN, index, 0, size()-1); + } #endif // Check for the existence of energies and energy widths - if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) + if (m_energies.size() != ebins() || m_ewidth.size() != ebins()) { throw GLATException::no_energies(G_SET_BIN); + } // Check for the existence of sky directions and solid angles - if (m_dirs.size() != npix() || m_omega.size() != npix()) + if (m_dirs.size() != npix() || m_omega.size() != npix()) { throw GLATException::no_dirs(G_SET_BIN); + } // Get pixel and energy bin indices. m_bin.m_index = index; @@ -904,10 +948,3 @@ void GLATEventCube::set_bin(const int& index) // Return return; } - - -/*========================================================================== - = = - = Friends = - = = - ==========================================================================*/ diff --git a/inst/lat/src/GLATException.cpp b/inst/lat/src/GLATException.cpp index 55b0e7c..f3f55ab 100644 --- a/inst/lat/src/GLATException.cpp +++ b/inst/lat/src/GLATException.cpp @@ -59,6 +59,24 @@ GLATException::file_open_error::file_open_error(std::string origin, /***********************************************************************//** + * @brief Member not set. + * + * @param[in] origin Name of method that has thrown the exception. + * @param[in] message Error message. + ***************************************************************************/ +GLATException::no_member::no_member(const std::string& origin, + const std::string& message) +{ + // Set origin and message + m_origin = origin; + m_message = message; + + // Return + return; +} + + +/***********************************************************************//** * @brief No sky pixels found. * * @param[in] origin Name of method that has thrown the exception. diff --git a/inst/lat/src/GLATResponseTable.cpp b/inst/lat/src/GLATResponseTable.cpp index 43864b1..96d8b76 100644 --- a/inst/lat/src/GLATResponseTable.cpp +++ b/inst/lat/src/GLATResponseTable.cpp @@ -21,7 +21,7 @@ /** * @file GLATResponseTable.cpp * @brief Fermi/LAT response table class implementation - * @author J. Knoedlseder + * @author Juergen Knoedlseder */ /* __ Includes ___________________________________________________________ */ @@ -159,7 +159,7 @@ void GLATResponseTable::clear(void) /***********************************************************************//** * @brief Clone instance -***************************************************************************/ + ***************************************************************************/ GLATResponseTable* GLATResponseTable::clone(void) const { return new GLATResponseTable(*this); @@ -290,9 +290,10 @@ int GLATResponseTable::index(const int& ie, const int& ic) const // Optionally check if the index is valid #if defined(G_RANGE_CHECK) if (ie < 0 || ie >= m_energy_num || - ic < 0 || ic >= m_ctheta_num) + ic < 0 || ic >= m_ctheta_num) { throw GException::out_of_range(G_INDEX, ie, ic, m_energy_num-1, m_ctheta_num-1); + } #endif // Compute index @@ -323,8 +324,9 @@ double GLATResponseTable::energy(const int& ie) const { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (ie < 0 || ie >= m_energy_num) + if (ie < 0 || ie >= m_energy_num) { throw GException::out_of_range(G_ENERGY, ie, 0, m_energy_num-1); + } #endif // Determine mean energy of bin @@ -462,8 +464,9 @@ double GLATResponseTable::energy_lo(const int& inx) const { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (inx < 0 || inx >= m_energy_num) + if (inx < 0 || inx >= m_energy_num) { throw GException::out_of_range(G_ENERGY_LO, inx, 0, m_energy_num-1); + } #endif // Return mean energy @@ -483,8 +486,9 @@ double GLATResponseTable::energy_hi(const int& inx) const { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (inx < 0 || inx >= m_energy_num) + if (inx < 0 || inx >= m_energy_num) { throw GException::out_of_range(G_ENERGY_HI, inx, 0, m_energy_num-1); + } #endif // Return mean energy @@ -504,8 +508,9 @@ double GLATResponseTable::costheta_lo(const int& inx) const { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (inx < 0 || inx >= m_ctheta_num) + if (inx < 0 || inx >= m_ctheta_num) { throw GException::out_of_range(G_COSTHETA_LO, inx, 0, m_ctheta_num-1); + } #endif // Return mean energy @@ -525,8 +530,9 @@ double GLATResponseTable::costheta_hi(const int& inx) const { // Optionally check if the index is valid #if defined(G_RANGE_CHECK) - if (inx < 0 || inx >= m_ctheta_num) + if (inx < 0 || inx >= m_ctheta_num) { throw GException::out_of_range(G_COSTHETA_HI, inx, 0, m_ctheta_num-1); + } #endif // Return mean energy @@ -562,10 +568,12 @@ std::vector GLATResponseTable::energies(void) const std::vector energies; // Push energies on vector - energies.push_back(m_energy[m_logE.inx_left()]); - energies.push_back(m_energy[m_logE.inx_left()]); - energies.push_back(m_energy[m_logE.inx_right()]); - energies.push_back(m_energy[m_logE.inx_right()]); + if (m_energy_num > 0) { + energies.push_back(m_energy[m_logE.inx_left()]); + energies.push_back(m_energy[m_logE.inx_left()]); + energies.push_back(m_energy[m_logE.inx_right()]); + energies.push_back(m_energy[m_logE.inx_right()]); + } // Return vector return energies; diff --git a/pyext/GEbounds.i b/pyext/GEbounds.i index 752d865..f73a2e7 100644 --- a/pyext/GEbounds.i +++ b/pyext/GEbounds.i @@ -66,7 +66,7 @@ public: GEnergy emax(int inx) const; GEnergy emean(int inx) const; GEnergy elogmean(int inx) const; - bool isin(const GEnergy& eng) const; + bool contains(const GEnergy& eng) const; }; diff --git a/src/fits/GFits.cpp b/src/fits/GFits.cpp index 9061889..4682ea7 100644 --- a/src/fits/GFits.cpp +++ b/src/fits/GFits.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * GFits.cpp - FITS file access class * + * GFits.cpp - FITS file access class * * ----------------------------------------------------------------------- * * copyright (C) 2008-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -348,10 +348,13 @@ void GFits::open(const std::string& filename, bool create) * @exception GException::fits_file_exist * Attemting to overwrite an existing file without having specified * clobber=true. + * @exception GException::fits_file_not_open + * FITS file needs to be opened before saving. + * + * Saves all HDUs of an open FITS file to disk. After saving, the FITS file + * remains open. Invoke the close() method if explicit closing is needed. + * Note that de-allocation of the GFits object also closes the FITS file. * - * Saves all HDUs to the FITS file by looping over the HDU's save() method. - * Saving does not mean that the FITS file is closed. Call the close() - * method if explicit closing is needed (de-allocation also closes the file). * In the special case that no first HDU exists an empty primary image is * created. ***************************************************************************/ @@ -368,6 +371,12 @@ void GFits::save(bool clobber) throw GException::fits_file_exist(G_SAVE, m_filename); } + // If no FITS file has been opened then throw an error + if (m_fitsfile == NULL) { + throw GException::fits_file_not_open(G_SAVE, + "Either open FITS file before saving or use saveto() method."); + } + // If no HDUs exist then save an empty primary image. if (size() == 0) { int status = 0; @@ -537,28 +546,33 @@ void GFits::append(const GFitsHDU& hdu) // Clone FITS HDU GFitsHDU* ptr = hdu.clone(); - // If FITS file exists then connect cloned HDU to FITS file - if (m_fitsfile != NULL) { - __fitsfile fptr = *(FPTR(m_fitsfile)); - fptr.HDUposition = n_hdu; - ptr->connect(&fptr); - } - else { - ptr->extno(n_hdu); - } + // Append HDU if it's valid + if (ptr != NULL) { - // Debug information - #if DEBUG - std::cout << "Append HDU (extno=" << n_hdu << ")." << std::endl; - #endif + // If FITS file exists then connect cloned HDU to FITS file + if (m_fitsfile != NULL) { + __fitsfile fptr = *(FPTR(m_fitsfile)); + fptr.HDUposition = n_hdu; + ptr->connect(&fptr); + } + else { + ptr->extno(n_hdu); + } - // Push back HDU - m_hdu.push_back(ptr); + // Debug information + #if DEBUG + std::cout << "Append HDU (extno=" << n_hdu << ")." << std::endl; + #endif - // Debug trailer - #if DEBUG - std::cout << "<-- GFits::append" << std::endl; - #endif + // Push back HDU + m_hdu.push_back(ptr); + + // Debug trailer + #if DEBUG + std::cout << "<-- GFits::append" << std::endl; + #endif + + } // endif: HDU was valid // Return return; diff --git a/src/fits/GFitsHDU.cpp b/src/fits/GFitsHDU.cpp index d7c7b6c..aea2102 100644 --- a/src/fits/GFitsHDU.cpp +++ b/src/fits/GFitsHDU.cpp @@ -358,6 +358,8 @@ void GFitsHDU::connect(void* vptr) /***********************************************************************//** * @brief Move FITS file pointer to HDU * + * @exception GException::fits_file_not_open + * No FITS file has been opened. * @exception GException::fits_hdu_not_found * Requested HDU not found. * @@ -366,6 +368,12 @@ void GFitsHDU::connect(void* vptr) ***************************************************************************/ void GFitsHDU::move_to_hdu(void) { + // Throw an exception if FITS file is not open + if (FPTR(m_fitsfile)->Fptr == NULL) { + throw GException::fits_file_not_open(G_MOVE_TO_HDU, + "Open file before moving file pointer to HDU."); + } + // Move to HDU int status = 0; status = __ffmahd(FPTR(m_fitsfile), m_hdunum+1, NULL, &status); @@ -382,13 +390,23 @@ void GFitsHDU::move_to_hdu(void) /***********************************************************************//** * @brief Get HDU type from FITS file * + * @exception GException::fits_file_not_open + * No FITS file has been opened. * @exception GException::fits_hdu_not_found * Requested HDU not found. ***************************************************************************/ GFitsHDU::HDUType GFitsHDU::get_hdu_type(void) const { + // Initialise HDU type + int type = 0; + + // Throw an exception if FITS file is not open + if (FPTR(m_fitsfile)->Fptr == NULL) { + throw GException::fits_file_not_open(G_GET_HDU_TYPE, + "Open file before requesting HDU type."); + } + // Get HDU type - int type = 0; int status = 0; status = __ffghdt(FPTR(m_fitsfile), &type, &status); if (status != 0) { diff --git a/src/fits/GFitsImage.cpp b/src/fits/GFitsImage.cpp index d9381de..b13e36b 100644 --- a/src/fits/GFitsImage.cpp +++ b/src/fits/GFitsImage.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * GFitsImage.cpp - FITS image class * + * GFitsImage.cpp - FITS image class * * ----------------------------------------------------------------------- * * copyright (C) 2008-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -371,8 +371,9 @@ int GFitsImage::naxes(int axis) const { // Check if axis is within the range #if defined(G_RANGE_CHECK) - if (axis < 0 || axis >= m_naxis) + if (axis < 0 || axis >= m_naxis) { throw GException::out_of_range(G_NAXES, axis, 0, m_naxis-1); + } #endif // Get axis dimension @@ -448,9 +449,10 @@ std::string GFitsImage::print(void) const result.append(parformat("Image type")+typecode(type())+"\n"); result.append(parformat("Number of dimensions")+str(naxis())+"\n"); result.append(parformat("Number of image pixels")+str(size())); - for (int i = 0; i < naxis(); ++i) + for (int i = 0; i < naxis(); ++i) { result.append("\n"+parformat("Number of bins in "+str(i)) + str(naxes(i))); + } // Append header information result.append(+"\n"+m_header.print()); @@ -577,6 +579,9 @@ void GFitsImage::init_image_header(void) * * @param[in] vptr FITS file void pointer. * + * @exception GException::fits_error + * FITS error. + * * Open FITS image in FITS file. Opening means connecting the FITS file * pointer to the image and reading the image and axes dimensions. ***************************************************************************/ @@ -585,8 +590,9 @@ void GFitsImage::open_image(void* vptr) // Move to HDU int status = 0; status = __ffmahd(FPTR(vptr), (FPTR(vptr)->HDUposition)+1, NULL, &status); - if (status != 0) + if (status != 0) { throw GException::fits_error(G_OPEN_IMAGE, status); + } // Save the FITS file pointer and the HDU number FPTR_COPY(m_fitsfile, vptr); @@ -594,8 +600,9 @@ void GFitsImage::open_image(void* vptr) // Get the image dimensions status = __ffgidm(FPTR(m_fitsfile), &m_naxis, &status); - if (status != 0) + if (status != 0) { throw GException::fits_error(G_OPEN_IMAGE, status); + } // Reset number of image pixels m_num_pixels = 0; @@ -609,13 +616,15 @@ void GFitsImage::open_image(void* vptr) // Get the axes dimensions status = __ffgisz(FPTR(m_fitsfile), m_naxis, m_naxes, &status); - if (status != 0) + if (status != 0) { throw GException::fits_error(G_OPEN_IMAGE, status); + } // Calculate number of image pixels m_num_pixels = 1; - for (int i = 0; i < m_naxis; ++i) + for (int i = 0; i < m_naxis; ++i) { m_num_pixels *= m_naxes[i]; + } } // endif: there is an image @@ -632,6 +641,9 @@ void GFitsImage::open_image(void* vptr) * @param[in] nulval Pointer to pixel nul value. * @param[out] anynul Number of nul values encountered during loading. * + * @exception GException::fits_error + * FITS error. + * * Load image pixels from FITS file. ***************************************************************************/ void GFitsImage::load_image(int datatype, const void* pixels, @@ -656,8 +668,9 @@ void GFitsImage::load_image(int datatype, const void* pixels, delete [] fpixel; delete [] lpixel; delete [] inc; - if (status != 0) + if (status != 0) { throw GException::fits_error(G_LOAD_IMAGE, status); + } } // Return @@ -671,12 +684,23 @@ void GFitsImage::load_image(int datatype, const void* pixels, * @param[in] datatype Datatype of pixels to be saved * @param[in] pixels Pixel array to be saved * + * @exception GException::fits_file_not_open + * No FITS file has been opened. + * @exception GException::fits_error + * FITS error. + * * Save image pixels into FITS file. In case that the HDU does not exist it * is created. In case that the pixel array is empty no data are saved; all * image pixels will be empty in this case. ***************************************************************************/ void GFitsImage::save_image(int datatype, const void* pixels) { + // Throw an exception if FITS file is not open + if (FPTR(m_fitsfile)->Fptr == NULL) { + throw GException::fits_file_not_open(G_SAVE_IMAGE, + "Open file before saving the image."); + } + // Move to HDU. We use here an explicit cfitsio moveto function since we // want to recover the error code ... int status = 0; @@ -686,23 +710,27 @@ void GFitsImage::save_image(int datatype, const void* pixels) if (status == 107) { status = 0; status = __ffcrim(FPTR(m_fitsfile), m_bitpix, m_naxis, m_naxes, &status); - if (status != 0) + if (status != 0) { throw GException::fits_error(G_SAVE_IMAGE, status); + } } - else if (status != 0) + else if (status != 0) { throw GException::fits_error(G_SAVE_IMAGE, status); + } // If HDU seems to be empty then create it now. This is only needed for the // primary HDU, since __ffmahd gives no error if the primary HDU is empty. // By checking the number of keywords in the HDU we detect an empty HDU ... int num = 0; status = __ffghsp(FPTR(m_fitsfile), &num, NULL, &status); - if (status != 0) + if (status != 0) { throw GException::fits_error(G_SAVE_IMAGE, status); + } if (num == 0) { status = __ffcrim(FPTR(m_fitsfile), m_bitpix, m_naxis, m_naxes, &status); - if (status != 0) + if (status != 0) { throw GException::fits_error(G_SAVE_IMAGE, status); + } } // Save the image pixels (if there are some ...) @@ -717,8 +745,9 @@ void GFitsImage::save_image(int datatype, const void* pixels) (void*)pixels, &status); delete [] fpixel; delete [] lpixel; - if (status != 0) + if (status != 0) { throw GException::fits_error(G_SAVE_IMAGE, status); + } } // Return @@ -749,8 +778,9 @@ void GFitsImage::fetch_data(void) init_data(); // If a FITS file is attached then load pixels from FITS file. - if (FPTR(m_fitsfile)->Fptr != NULL) + if (FPTR(m_fitsfile)->Fptr != NULL) { load_image(type(), ptr_data(), ptr_nulval(), &m_anynul); + } } // endif: there were pixels available @@ -771,8 +801,9 @@ int GFitsImage::offset(const int& ix) const { // Check if axis is within the range #if defined(G_RANGE_CHECK) - if (ix < 0 || ix >= m_num_pixels) + if (ix < 0 || ix >= m_num_pixels) { throw GException::out_of_range(G_OFFSET_1D, ix, 0, m_num_pixels-1); + } #endif // Return index @@ -797,15 +828,18 @@ int GFitsImage::offset(const int& ix) const int GFitsImage::offset(const int& ix, const int& iy) const { // Operator is only valid for 2D images - if (m_naxis < 2) + if (m_naxis < 2) { throw GException::fits_wrong_image_operator(G_OFFSET_2D, m_naxis, 2); + } // Check if axis is within the range #if defined(G_RANGE_CHECK) - if (ix < 0 || ix >= m_naxes[0]) + if (ix < 0 || ix >= m_naxes[0]) { throw GException::out_of_range(G_OFFSET_2D, ix, 0, m_naxes[0]-1); - if (iy < 0 || iy >= m_naxes[1]) + } + if (iy < 0 || iy >= m_naxes[1]) { throw GException::out_of_range(G_OFFSET_2D, iy, 0, m_naxes[1]-1); + } #endif // Return offset @@ -831,17 +865,21 @@ int GFitsImage::offset(const int& ix, const int& iy) const int GFitsImage::offset(const int& ix, const int& iy, const int& iz) const { // Operator is only valid for 3D images - if (m_naxis < 3) + if (m_naxis < 3) { throw GException::fits_wrong_image_operator(G_OFFSET_3D, m_naxis, 3); + } // Check if axis is within the range #if defined(G_RANGE_CHECK) - if (ix < 0 || ix >= m_naxes[0]) + if (ix < 0 || ix >= m_naxes[0]) { throw GException::out_of_range(G_OFFSET_3D, ix, 0, m_naxes[0]-1); - if (iy < 0 || iy >= m_naxes[1]) + } + if (iy < 0 || iy >= m_naxes[1]) { throw GException::out_of_range(G_OFFSET_3D, iy, 0, m_naxes[1]-1); - if (iz < 0 || iz >= m_naxes[2]) + } + if (iz < 0 || iz >= m_naxes[2]) { throw GException::out_of_range(G_OFFSET_3D, iz, 0, m_naxes[2]-1); + } #endif // Return offset @@ -869,19 +907,24 @@ int GFitsImage::offset(const int& ix, const int& iy, const int& iz, const int& it) const { // Operator is only valid for 4D images - if (m_naxis < 4) + if (m_naxis < 4) { throw GException::fits_wrong_image_operator(G_OFFSET_4D, m_naxis, 4); + } // Check if axis is within the range #if defined(G_RANGE_CHECK) - if (ix < 0 || ix >= m_naxes[0]) + if (ix < 0 || ix >= m_naxes[0]) { throw GException::out_of_range(G_OFFSET_4D, ix, 0, m_naxes[0]-1); - if (iy < 0 || iy >= m_naxes[1]) + } + if (iy < 0 || iy >= m_naxes[1]) { throw GException::out_of_range(G_OFFSET_4D, iy, 0, m_naxes[1]-1); - if (iz < 0 || iz >= m_naxes[2]) + } + if (iz < 0 || iz >= m_naxes[2]) { throw GException::out_of_range(G_OFFSET_4D, iz, 0, m_naxes[2]-1); - if (it < 0 || it >= m_naxes[3]) + } + if (it < 0 || it >= m_naxes[3]) { throw GException::out_of_range(G_OFFSET_4D, it, 0, m_naxes[3]-1); + } #endif // Return offset @@ -929,8 +972,9 @@ void GFitsImage::copy_members(const GFitsImage& image) m_naxes = NULL; if (image.m_naxes != NULL && m_naxis > 0) { m_naxes = new long[m_naxis]; - for (int i = 0; i < m_naxis; ++i) + for (int i = 0; i < m_naxis; ++i) { m_naxes[i] = image.m_naxes[i]; + } } // Return diff --git a/src/model/GModels.cpp b/src/model/GModels.cpp index acd8939..0d26de9 100644 --- a/src/model/GModels.cpp +++ b/src/model/GModels.cpp @@ -312,8 +312,13 @@ GModels* GModels::clone(void) const ***************************************************************************/ void GModels::append(const GModel& model) { - // Append model - m_models.push_back(model.clone()); + // Make deep copy of model + GModel* model_copy = model.clone(); + + // If copy is valid, append it to the container + if (model_copy != NULL) { + m_models.push_back(model_copy); + } // Set parameter pointers set_pointers(); diff --git a/src/obs/GEbounds.cpp b/src/obs/GEbounds.cpp index 7b88264..0530a8a 100644 --- a/src/obs/GEbounds.cpp +++ b/src/obs/GEbounds.cpp @@ -578,7 +578,7 @@ GEnergy GEbounds::elogmean(int inx) const * * @param[in] eng Energy to be tested. ***************************************************************************/ -bool GEbounds::isin(const GEnergy& eng) const +bool GEbounds::contains(const GEnergy& eng) const { // Initialise test bool found = false; diff --git a/src/obs/GEnergy.cpp b/src/obs/GEnergy.cpp index f28b401..2024ef1 100644 --- a/src/obs/GEnergy.cpp +++ b/src/obs/GEnergy.cpp @@ -228,7 +228,7 @@ double GEnergy::TeV(void) const double GEnergy::log10keV(void) const { // Return log10 energy - return (log10keV()+3.0); + return (log10MeV()+3.0); } diff --git a/src/opt/GOptimizerLM.cpp b/src/opt/GOptimizerLM.cpp index 3900297..53357f0 100644 --- a/src/opt/GOptimizerLM.cpp +++ b/src/opt/GOptimizerLM.cpp @@ -1,7 +1,7 @@ /*************************************************************************** - * GOptimizerLM.cpp - Levenberg Marquardt optimizer * + * GOptimizerLM.cpp - Levenberg Marquardt optimizer * * ----------------------------------------------------------------------- * - * copyright (C) 2009-2011 by Jurgen Knodlseder * + * copyright (C) 2009-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file GOptimizerLM.cpp * @brief Levenberg-Marquardt optimizer class implementation - * @author J. Knodlseder + * @author Juergen Knoedlseder */ /* __ Includes ___________________________________________________________ */ @@ -948,108 +948,108 @@ double GOptimizerLM::step_size(GVector* grad, GOptimizerPars* pars) ***************************************************************************/ void GOptimizerLM::errors(GOptimizerFunction* fct, GOptimizerPars* pars) { - // Get number of parameters - int npars = pars->npars(); + // Continue only if pointers are valid + if (fct != NULL && pars != NULL) { - // Perform final parameter evaluation - fct->eval(*pars); + // Get number of parameters + int npars = pars->npars(); - // Fetch sparse matrix pointer. We have to do this after the eval() - // method since eval() will allocate new memory for the covariance - // matrix! - GSparseMatrix* covar = fct->covar(); + // Perform final parameter evaluation + fct->eval(*pars); - // Save best fitting value - m_value = *(fct->value()); + // Fetch sparse matrix pointer. We have to do this after the eval() + // method since eval() will allocate new memory for the covariance + // matrix! + GSparseMatrix* covar = fct->covar(); - // Save covariance matrix - GSparseMatrix save_covar = GSparseMatrix(*covar); + // Save best fitting value + m_value = *(fct->value()); - // Signal no diagonal element loading - bool diag_loaded = false; + // Save covariance matrix + GSparseMatrix save_covar = GSparseMatrix(*covar); - // Loop over error computation (maximum 2 turns) - for (int i = 0; i < 2; ++i) { + // Signal no diagonal element loading + bool diag_loaded = false; - // Solve: covar * X = unit - try { - covar->cholesky_decompose(1); - GVector unit(npars); - for (int ipar = 0; ipar < npars; ++ipar) { - unit[ipar] = 1.0; - GVector x = covar->cholesky_solver(unit,1); - if (x[ipar] >= 0.0) - pars->par(ipar).error(sqrt(x[ipar])); - else { - pars->par(ipar).error(0.0); - m_status = G_LM_BAD_ERRORS; + // Loop over error computation (maximum 2 turns) + for (int i = 0; i < 2; ++i) { + + // Solve: covar * X = unit + try { + covar->cholesky_decompose(1); + GVector unit(npars); + for (int ipar = 0; ipar < npars; ++ipar) { + unit[ipar] = 1.0; + GVector x = covar->cholesky_solver(unit,1); + if (x[ipar] >= 0.0) { + pars->par(ipar).error(sqrt(x[ipar])); + } + else { + pars->par(ipar).error(0.0); + m_status = G_LM_BAD_ERRORS; + } + unit[ipar] = 0.0; } - unit[ipar] = 0.0; - } - } - catch (GException::matrix_zero &e) { - m_status = G_LM_SINGULAR; - if (m_logger != NULL) { - *m_logger << "GOptimizerLM::terminate: " - << "All curvature matrix elements are zero." - << std::endl; } - break; - } - catch (GException::matrix_not_pos_definite &e) { - - // Load diagonal if this has not yet been tried - if (!diag_loaded) { - - // Flag errors as inaccurate - m_status = G_LM_BAD_ERRORS; + catch (GException::matrix_zero &e) { + m_status = G_LM_SINGULAR; if (m_logger != NULL) { - *m_logger << "Non-Positive definite curvature matrix encountered." - << std::endl; - *m_logger << "Load diagonal elements with 1e-10." - << " Fit errors may be inaccurate." + *m_logger << "GOptimizerLM::terminate: " + << "All curvature matrix elements are zero." << std::endl; } + break; + } + catch (GException::matrix_not_pos_definite &e) { - // Try now with diagonal loaded matrix - *covar = save_covar; - for (int ipar = 0; ipar < npars; ++ipar) - (*covar)(ipar,ipar) += 1.0e-10; + // Load diagonal if this has not yet been tried + if (!diag_loaded) { - // Signal loading - diag_loaded = true; + // Flag errors as inaccurate + m_status = G_LM_BAD_ERRORS; + if (m_logger != NULL) { + *m_logger << "Non-Positive definite curvature matrix encountered." + << std::endl; + *m_logger << "Load diagonal elements with 1e-10." + << " Fit errors may be inaccurate." + << std::endl; + } - // Try again - continue; + // Try now with diagonal loaded matrix + *covar = save_covar; + for (int ipar = 0; ipar < npars; ++ipar) { + (*covar)(ipar,ipar) += 1.0e-10; + } - } // endif: diagonal has not yet been loaded + // Signal loading + diag_loaded = true; - // ... otherwise signal an error - else { - m_status = G_LM_NOT_POSTIVE_DEFINITE; - if (m_logger != NULL) { - *m_logger << "Non-Positive definite curvature matrix encountered," - << " even after diagonal loading." << std::endl; + // Try again + continue; + + } // endif: diagonal has not yet been loaded + + // ... otherwise signal an error + else { + m_status = G_LM_NOT_POSTIVE_DEFINITE; + if (m_logger != NULL) { + *m_logger << "Non-Positive definite curvature matrix encountered," + << " even after diagonal loading." << std::endl; + } + break; } - break; } - } - catch (std::exception &e) { - throw; - } + catch (std::exception &e) { + throw; + } - // If no error occured then break now - break; + // If no error occured then break now + break; + + } // endfor: looped over error computation - } // endfor: looped over error computation + } // endif: pointers were valid // Return return; } - - -/*========================================================================== - = = - = Friends = - = = - ==========================================================================*/ diff --git a/src/sky/GException_sky.cpp b/src/sky/GException_sky.cpp index 7b47989..74e4f51 100644 --- a/src/sky/GException_sky.cpp +++ b/src/sky/GException_sky.cpp @@ -255,7 +255,7 @@ GException::wcs_hpx_bad_ordering::wcs_hpx_bad_ordering(std::string origin, * @brief Singular matrix encountered * * @param[in] origin Method that throws the error. - * @param[in] naxis Number of matrix axes. + * @param[in] naxis Matrix size (naxis*naxis). * @param[in] mat Matrix (stored in vector). ***************************************************************************/ GException::wcs_singular_matrix::wcs_singular_matrix(std::string origin, @@ -266,13 +266,20 @@ GException::wcs_singular_matrix::wcs_singular_matrix(std::string m_origin = origin; m_message = "Singular matrix encountered: mat=["; for (int i = 0, ij = 0; i < naxis; ++i) { - if (i > 0) + if (i > 0) { m_message += ","; + } m_message += "["; for (int j = 0; j < naxis; ++j, ++ij) { - if (j > 0) + if (j > 0) { m_message += ","; - m_message += str(mat[ij]); + } + if (ij < mat.size()) { + m_message += str(mat[ij]); + } + else { + m_message += " ... "; + } } m_message += "]"; } diff --git a/src/support/GNodeArray.cpp b/src/support/GNodeArray.cpp index 819f706..da44b0f 100644 --- a/src/support/GNodeArray.cpp +++ b/src/support/GNodeArray.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * GNodeArray.cpp - Array of nodes class * + * GNodeArray.cpp - Array of nodes class * * ----------------------------------------------------------------------- * * copyright (C) 2008-2012 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * @@ -35,7 +35,9 @@ /* __ Method name definitions ____________________________________________ */ #define G_ACCESS "GNodeArray::operator[](int)" -#define G_INTERPOLATE "GNodeArray::interpolate(std::vector&, double&" +#define G_INTERPOLATE "GNodeArray::interpolate(double&,std::vector&)" +#define G_SET_VALUE "GNodeArray::set_value(double&)" +#define G_SETUP "GNodeArray::setup()" /* __ Macros _____________________________________________________________ */ @@ -166,7 +168,8 @@ GNodeArray::~GNodeArray(void) /***********************************************************************//** * @brief Assignment operator * - * @param[in] array Instance that should be assigned. + * @param[in] array Node array. + * @return Node array. ***************************************************************************/ GNodeArray& GNodeArray::operator= (const GNodeArray& array) { @@ -193,6 +196,7 @@ GNodeArray& GNodeArray::operator= (const GNodeArray& array) * @brief Node access operator * * @param[in] index Node index [0,...,size()-1]. + * @return Node value. * * @exception GException::out_of_range * Node index is out of range. @@ -201,8 +205,9 @@ double& GNodeArray::operator[](const int& index) { // Compile option: raise an exception if index is out of range #if defined(G_RANGE_CHECK) - if (index < 0 || index >= size()) + if (index < 0 || index >= size()) { throw GException::out_of_range(G_ACCESS, index, size()-1); + } #endif // Return node @@ -214,6 +219,7 @@ double& GNodeArray::operator[](const int& index) * @brief Node access operator (const version) * * @param[in] index Node index [0,...,size()-1]. + * @return Node value. * * @exception GException::out_of_range * Node index is out of range. @@ -222,8 +228,9 @@ const double& GNodeArray::operator[](const int& index) const { // Compile option: raise an exception if index is out of range #if defined(G_RANGE_CHECK) - if (index < 0 || index >= size()) + if (index < 0 || index >= size()) { throw GException::out_of_range(G_ACCESS, index, size()-1); + } #endif // Return node @@ -239,7 +246,7 @@ const double& GNodeArray::operator[](const int& index) const /***********************************************************************//** * @brief Clear instance -***************************************************************************/ + ***************************************************************************/ void GNodeArray::clear(void) { // Free class members @@ -255,7 +262,9 @@ void GNodeArray::clear(void) /***********************************************************************//** * @brief Clone instance -***************************************************************************/ + * + * @return Deep copy of node array. + ***************************************************************************/ GNodeArray* GNodeArray::clone(void) const { return new GNodeArray(*this); @@ -267,6 +276,8 @@ GNodeArray* GNodeArray::clone(void) const * * @param[in] num Number of nodes * @param[in] array Node values \f$x_i\f$. + * + * Setup node array from an array of double precision values. ***************************************************************************/ void GNodeArray::nodes(const int& num, const double* array) { @@ -277,8 +288,9 @@ void GNodeArray::nodes(const int& num, const double* array) init_members(); // Set node values - for (int i = 0; i < num; ++i) + for (int i = 0; i < num; ++i) { m_node.push_back(array[i]); + } // Setup node distances and linear array handling setup(); @@ -292,6 +304,8 @@ void GNodeArray::nodes(const int& num, const double* array) * @brief Set node array from vector * * @param[in] vector Vector from which node array will be built. + * + * Setup node array from a vector of values. ***************************************************************************/ void GNodeArray::nodes(const GVector& vector) { @@ -302,8 +316,9 @@ void GNodeArray::nodes(const GVector& vector) init_members(); // Set node values - for (int i = 0; i < vector.size(); ++i) + for (int i = 0; i < vector.size(); ++i) { m_node.push_back(vector[i]); + } // Setup node distances and linear array handling setup(); @@ -317,6 +332,8 @@ void GNodeArray::nodes(const GVector& vector) * @brief Set node array from vector * * @param[in] vector Vector from which node array will be built. + * + * Setup node array from a vector of double precision values. ***************************************************************************/ void GNodeArray::nodes(const std::vector& vector) { @@ -373,13 +390,15 @@ double GNodeArray::interpolate(const double& value, const std::vector& vector) const { // Throw exception if there are not enough nodes - if (m_node.size() < 2) + if (m_node.size() < 2) { throw GException::not_enough_nodes(G_INTERPOLATE, m_node.size()); + } // Throw exception if vectors have not the same size - if (m_node.size() != vector.size()) + if (m_node.size() != vector.size()) { throw GException::vector_mismatch(G_INTERPOLATE, m_node.size(), vector.size()); + } // Set interpolation value (circumvent const correctness) const_cast(this)->set_value(value); @@ -398,6 +417,9 @@ double GNodeArray::interpolate(const double& value, * * @param[in] value Value for which the interpolation should be done. * + * @exception GException::not_enough_nodes + * At least two nodes are required for setting up the factors + * * Set the indices that bound the specified value and the corresponding * weighting factors for linear interpolation. If the array has a linear * form (i.e. the nodes are equidistant), an analytic formula is used to @@ -408,55 +430,60 @@ void GNodeArray::set_value(const double& value) { // Get number of nodes int nodes = m_node.size(); + + // Throw an exception if less than 2 nodes are available + if (nodes < 2) { + throw GException::not_enough_nodes(G_SET_VALUE, nodes); + } - // Continue only if we have at least 2 nodes - if (nodes > 1) { - - // If array is linear then get left index from analytic formula - if (m_is_linear) { + // If array is linear then get left index from analytic formula + if (m_is_linear) { - // Set left index - m_inx_left = int(m_linear_slope * value + m_linear_offset); + // Set left index + m_inx_left = int(m_linear_slope * value + m_linear_offset); + + // Keep index in valid range + if (m_inx_left < 0) m_inx_left = 0; + else if (m_inx_left >= nodes-1) m_inx_left = nodes - 2; - // Keep index in valid range - if (m_inx_left < 0) m_inx_left = 0; - else if (m_inx_left >= nodes-1) m_inx_left = nodes - 2; + } // endif: array is linear - } // endif: array is linear + // ... otherwise search the relevant indices by bisection + else { + + // Set left index if value is before first node + if (value < m_node[0]) { + m_inx_left = 0; + } - // ... otherwise search the relevant indices by bisection + // Set left index if value is after last node + else if (value > m_node[nodes-1]) { + m_inx_left = nodes - 2; + } + + // Set left index by bisection else { - // Set left index if value is before first node - if (value < m_node[0]) - m_inx_left = 0; - - // Set left index if value is after last node - else if (value > m_node[nodes-1]) - m_inx_left = nodes - 2; - - // Set left index by bisection - else { - int low = 0; - int high = nodes - 1; - while ((high - low) > 1) { - int mid = (low+high) / 2; - if (m_node[mid] > value) - high = mid; - else if (m_node[mid] <= value) - low = mid; + int low = 0; + int high = nodes - 1; + while ((high - low) > 1) { + int mid = (low+high) / 2; + if (m_node[mid] > value) { + high = mid; } - m_inx_left = low; - } // endelse: did bisection - } + else if (m_node[mid] <= value) { + low = mid; + } + } + m_inx_left = low; + } // endelse: did bisection + } - // Set right index - m_inx_right = m_inx_left + 1; + // Set right index + m_inx_right = m_inx_left + 1; - // Set weighting factors - m_wgt_right = (value - m_node[m_inx_left]) / m_step[m_inx_left]; - m_wgt_left = 1.0 - m_wgt_right; - - } // endif: there were at least 2 nodes + // Set weighting factors + m_wgt_right = (value - m_node[m_inx_left]) / m_step[m_inx_left]; + m_wgt_left = 1.0 - m_wgt_right; // Return return; @@ -493,7 +520,9 @@ std::string GNodeArray::print(void) const for (int i = 0; i < size(); ++i) { result.append("\n"+parformat("Node "+str(i))); result.append(str(m_node[i])); - result.append(" (delta="+str(m_step[i])+")"); + if (i < m_step.size()) { + result.append(" (delta="+str(m_step[i])+")"); + } } // Return result @@ -515,7 +544,7 @@ void GNodeArray::init_members(void) // Initialise members m_node.clear(); m_step.clear(); - m_is_linear = 0; + m_is_linear = false; m_linear_slope = 0.0; m_linear_offset = 0.0; m_inx_left = 0; @@ -563,6 +592,11 @@ void GNodeArray::free_members(void) /***********************************************************************//** * @brief Compute distance array and linear slope/offset + * + * Precomputes values for fast interpolation. The precomputation requires + * at least 2 nodes to be present in the node array. If less then two + * nodes are present, the distance vector m_step will be empty and no + * computation is done. ***************************************************************************/ void GNodeArray::setup(void) { @@ -571,13 +605,14 @@ void GNodeArray::setup(void) // Get number of nodes int nodes = m_node.size(); - - // Continue only if we have nodes at least 2 nodes + + // Continue only if there are at least 2 nodes if (nodes > 1) { - + // Setup distance array between subsequent nodes - for (int i = 0; i < nodes-1; ++i) + for (int i = 0; i < nodes-1; ++i) { m_step.push_back(m_node[i+1] - m_node[i]); + } // Evaluate linear slope and offset m_linear_slope = double(nodes-1) / (m_node[nodes-1] - m_node[0]); @@ -592,9 +627,10 @@ void GNodeArray::setup(void) break; } } - - } // endif: there were at least 2 nodes + + } // endif: there were at least two nodes // Return return; } + diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index 5aadebb..405453b 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -258,8 +258,8 @@ void TestGSky::test_GWcslib(void) try { // Allocate projection CEL and GAL - GWcslib *cel = (GWcslib*)registry.alloc(registry.code(i)); - GWcslib *gal = (GWcslib*)registry.alloc(registry.code(i)); + GWcslib *cel = dynamic_cast(registry.alloc(registry.code(i))); + GWcslib *gal = dynamic_cast(registry.alloc(registry.code(i))); // Throw an error if allocation failed if (cel == NULL || gal == NULL) {