diff --git a/.gitignore b/.gitignore index 7f33671..32d1c86 100644 --- a/.gitignore +++ b/.gitignore @@ -54,7 +54,7 @@ aclocal.m4 autom4te.cache config.guess config.h -config.h.in +config.h.in* config.status config.sub configure @@ -65,6 +65,7 @@ test/test_url.dat test/*.xml test/reports/*.xml test/*.trs +test/region.ds9 examples/python/*.xml examples/python/*.html @@ -82,9 +83,11 @@ test/test_GSky test/test_GSupport test/test_GVector test/test_GXml +test/test_GXspec test/test_GVO test/test_LAT test/test_MWL test/test_CTA test/test_COM examples/cpp/readmodel/readmodel + diff --git a/AUTHORS b/AUTHORS index e461a2c..773fcde 100644 --- a/AUTHORS +++ b/AUTHORS @@ -1,4 +1,6 @@ Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu) Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr) Christoph Deil (Christoph.Deil@mpi-hd.mpg.de) -Michael Mayer (michael.mayer@desy.de) \ No newline at end of file +Michael Mayer (michael.mayer@desy.de) +Anneli Schulz (anneli.schulz@desy.de) +Chia-Chun Lu (Chia-Chun.Lu@mpi-hd.mpg.de) \ No newline at end of file diff --git a/ChangeLog b/ChangeLog index 9b891c1..9f9a8e1 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2013-09-07 Juergen Knoedlseder +2013-11-01 Juergen Knoedlseder * Version 0.8.0 released ======================== @@ -12,6 +12,9 @@ Optimize computations Introduce gammalib namespace for constants and functions GModelSpatialDiffuseCube reads ENERGIES extension + Add sky region handling + Add on/off analysis support for CTA + Add xspec interface support 2012-12-18 Juergen Knoedlseder diff --git a/NEWS b/NEWS index 2398049..1aa6202 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 00-08-00 Author(s) : Juergen Knoedlseder -September, 7th, 2013 +1 November 2013 1. Introduction @@ -221,6 +221,12 @@ Add a sin() and cos() pre computation cache to speed up repeated distance computations using dist() or dist_deg(). Coordinate access methods are now declared inline to provide also faster data access. +Classes have been added to support handling of sky regions. The +GSkyRegion class is the abstract base class of all sky regions. The +GSkyRegionCircle class implements a circle on the sky. The DS9 region +format is used for loading and saving region information into a file. +Regions a gathered in a regions container, which is implement by the +GSkyRegions class. 13. Support module ------------------ @@ -251,29 +257,41 @@ with a SAMP Hub. The method supports SAMP discovery and registering and unregistering. -16. Test module +16. Xspec module +---------------- +A module has been added to provide an interface to the Xspec analysis +package. Three classes have been implemented that support the Xspec +file formats: GPha, GArf and GRmf. + + +17. Test module --------------- None -17. Multi-wavelength interface +18. Multi-wavelength interface ------------------------------ None -18. Fermi-LAT interface +19. Fermi-LAT interface ---------------------- None -19. CTA interface +20. CTA interface ----------------- The computations for diffuse map templates have been optimized by introducing a pre computation cache for IRF and Npred values. This considerably reduces the computing time for diffuse spatial maps. +Support for on/off analysis has been added. The GCTAOnOffObservation +class gathers all information that is required for an on/off analysis +for a given run. The GCTAOnOffObservations container allows to analyse +several runs in one shot. + -20. COMPTEL interface +21. COMPTEL interface --------------------- None diff --git a/configure.ac b/configure.ac index 5042a9d..132f18a 100644 --- a/configure.ac +++ b/configure.ac @@ -860,7 +860,8 @@ else pyext/gammalib/sky_wrap.cpp pyext/gammalib/sky.py pyext/gammalib/support_wrap.cpp pyext/gammalib/support.py pyext/gammalib/test_wrap.cpp pyext/gammalib/test.py - pyext/gammalib/xml_wrap.cpp pyext/gammalib/xml.py],, + pyext/gammalib/xml_wrap.cpp pyext/gammalib/xml.py + pyext/gammalib/xspec_wrap.cpp pyext/gammalib/xspec.py],, [has_wrappers="no"]) if test "x$with_mwl" = "xyes"; then AC_CHECK_FILES([pyext/gammalib/mwl_wrap.cpp pyext/gammalib/mwl.py],, @@ -1239,6 +1240,7 @@ AC_CONFIG_FILES([Makefile src/fits/Makefile src/xml/Makefile src/vo/Makefile + src/xspec/Makefile src/sky/Makefile src/opt/Makefile src/obs/Makefile diff --git a/examples/cpp/readmodel/Makefile b/examples/cpp/readmodel/Makefile index fa7f9b5..99ec394 100644 --- a/examples/cpp/readmodel/Makefile +++ b/examples/cpp/readmodel/Makefile @@ -4,9 +4,5 @@ LDFLAGS=-L/usr/local/gamma/lib -lgamma DEPS= OBJ=readmodel.cpp -%o: *.cpp $(DEPS) - $(CXX) -c -o $@ $< $(CFLAGS) - readmodel: $(OBJ) $(CXX) -o $@ $^ $(CFLAGS) $(LDFLAGS) - diff --git a/include/GArf.hpp b/include/GArf.hpp new file mode 100644 index 0000000..9d30270 --- /dev/null +++ b/include/GArf.hpp @@ -0,0 +1,157 @@ +/*************************************************************************** + * GArf.hpp - XSPEC Auxiliary Response File class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GArf.hpp + * @brief XSPEC Auxiliary Response File class definition + * @author Juergen Knoedlseder + */ + +#ifndef GARF_HPP +#define GARF_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GBase.hpp" +#include "GEbounds.hpp" +#include "GFits.hpp" +#include "GFitsTable.hpp" + + +/***********************************************************************//** + * @class GArf + * + * @brief Auxiliary Response File class + ***************************************************************************/ +class GArf : public GBase { + +public: + // Constructors and destructors + GArf(void); + explicit GArf(const std::string& filename); + explicit GArf(const GEbounds& ebds); + GArf(const GArf& arf); + virtual ~GArf(void); + + // Operators + GArf& operator=(const GArf& arf); + double& operator[](const int& index); + const double& operator[](const int& index) const; + + // Methods + void clear(void); + GArf* clone(void) const; + int size(void) const; + double& at(const int& index); + 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 write(GFits& fits) const; + const std::string& filename(void) const; + std::string print(const GChatter& chatter = NORMAL) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GArf& pha); + void free_members(void); + + // Protected members + mutable std::string m_filename; //!< Filename of origin + GEbounds m_ebounds; //!< Energy boundaries + std::vector m_specresp; //!< Spectral response +}; + + +/***********************************************************************//** + * @brief Return content of spectral bin + * + * @param[in] index Bin index [0,...,size()-1]. + * + * Returns reference to content of spectral bin with specified @p index. + ***************************************************************************/ +inline +double& GArf::operator[](const int& index) +{ + return (m_specresp[index]); +} + + +/***********************************************************************//** + * @brief Return content of spectral bin (const version) + * + * @param[in] index Bin index [0,...,size()-1]. + * + * Returns reference to content of spectral bin with specified @p index. + ***************************************************************************/ +inline +const double& GArf::operator[](const int& index) const +{ + return (m_specresp[index]); +} + + +/***********************************************************************//** + * @brief Return number of spectral bins + * + * @return Number of spectral bins. + * + * Returns the number of spectral bins. + ***************************************************************************/ +inline +int GArf::size(void) const +{ + return (m_specresp.size()); +} + + +/***********************************************************************//** + * @brief Return energy boundaries + * + * @return Energy boundaries for all spectral bins. + * + * Returns the energy boundaries for all spectral bins. + ***************************************************************************/ +inline +const GEbounds& GArf::ebounds(void) const +{ + return m_ebounds; +} + + +/***********************************************************************//** + * @brief Return file name + * + * @return File name from which the ARF information has been read or into + * which ARF information has been saved. + * + * Returns the file name from which the ARF information has been read or into + * which ARF information has been saved. The returned string will be empty if + * no load() or save() method has been called before. + ***************************************************************************/ +inline +const std::string& GArf::filename(void) const +{ + return (m_filename); +} + +#endif /* GARF_HPP */ diff --git a/include/GPha.hpp b/include/GPha.hpp new file mode 100644 index 0000000..ab2d5a4 --- /dev/null +++ b/include/GPha.hpp @@ -0,0 +1,210 @@ +/*************************************************************************** + * GPha.hpp - XSPEC Pulse Height Analyzer class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GPha.hpp + * @brief XSPEC Pulse Height Analyzer class definition + * @author Juergen Knoedlseder + */ + +#ifndef GPHA_HPP +#define GPHA_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GBase.hpp" +#include "GEnergy.hpp" +#include "GEbounds.hpp" +#include "GFits.hpp" +#include "GFitsTable.hpp" + + +/***********************************************************************//** + * @class GPha + * + * @brief Pulse Height Analyzer class + ***************************************************************************/ +class GPha : public GBase { + +public: + // Constructors and destructors + GPha(void); + explicit GPha(const std::string& filename); + explicit GPha(const GEbounds& ebds); + GPha(const GPha& pha); + virtual ~GPha(void); + + // Operators + GPha& operator=(const GPha& pha); + double& operator[](const int& index); + const double& operator[](const int& index) const; + + // Methods + void clear(void); + GPha* clone(void) const; + int size(void) const; + double& at(const int& index); + const double& at(const int& index) const; + const GEbounds& ebounds(void) const; + double counts(void) const; + const double& underflow(void) const; + const double& overflow(void) const; + 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 write(GFits& fits) const; + const std::string& filename(void) const; + std::string print(const GChatter& chatter = NORMAL) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GPha& pha); + void free_members(void); + + // Protected members + mutable std::string m_filename; //!< Filename of origin + GEbounds m_ebounds; //!< Energy boundaries + std::vector m_counts; //!< Counts data + double m_underflow; //!< Number of underflowing events + double m_overflow; //!< Number of overflowing events + double m_outflow; //!< Number of outflowing events +}; + + +/***********************************************************************//** + * @brief Return content of spectral bin + * + * @param[in] index Bin index [0,...,size()-1]. + * + * Returns reference to content of spectral bin with specified @p index. + ***************************************************************************/ +inline +double& GPha::operator[](const int& index) +{ + return (m_counts[index]); +} + + +/***********************************************************************//** + * @brief Return content of spectral bin (const version) + * + * @param[in] index Bin index [0,...,size()-1]. + * + * Returns reference to content of spectral bin with specified @p index. + ***************************************************************************/ +inline +const double& GPha::operator[](const int& index) const +{ + return (m_counts[index]); +} + + +/***********************************************************************//** + * @brief Return number of bins in spectrum + * + * @return Number of bins in spectrum. + * + * Returns the number of bins in the spectrum. + ***************************************************************************/ +inline +int GPha::size(void) const +{ + return (m_counts.size()); +} + + +/***********************************************************************//** + * @brief Return energy boundaries + * + * @return Energy boundaries for all spectral bins. + * + * Returns the energy boundaries for all spectral bins. + ***************************************************************************/ +inline +const GEbounds& GPha::ebounds(void) const +{ + return m_ebounds; +} + + +/***********************************************************************//** + * @brief Return number of underflow counts + * + * @return Number of counts with energies lower than the first energy bin. + * + * Returns the number of counts with energies lower than the first energy + * bin. + ***************************************************************************/ +inline +const double& GPha::underflow(void) const +{ + return m_underflow; +} + + +/***********************************************************************//** + * @brief Return number of overflow counts + * + * @return Number of counts with energies larger than the last energy bin. + * + * Returns the number of counts with energies larger than the last energy + * bin. + ***************************************************************************/ +inline +const double& GPha::overflow(void) const +{ + return m_overflow; +} + + +/***********************************************************************//** + * @brief Return number of outflow counts + * + * @return Number of counts with energies between energy bins. + * + * Returns the number of counts with energies between energy bins. + ***************************************************************************/ +inline +const double& GPha::outflow(void) const +{ + return m_outflow; +} + + +/***********************************************************************//** + * @brief Return file name + * + * @return File name from which the PHA information has been read or into + * which PHA information has been saved. + * + * Returns the file name from which the PHA information has been read or into + * which PHA information has been saved. The returned string will be empty if + * no load() or save() method has been called before. + ***************************************************************************/ +inline +const std::string& GPha::filename(void) const +{ + return (m_filename); +} + +#endif /* GPHA_HPP */ diff --git a/include/GRmf.hpp b/include/GRmf.hpp new file mode 100644 index 0000000..6f9b060 --- /dev/null +++ b/include/GRmf.hpp @@ -0,0 +1,211 @@ +/*************************************************************************** + * GRmf.hpp - XSPEC Redistribution Matrix File class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GRmf.hpp + * @brief XSPEC Redistribution Matrix File class definition + * @author Juergen Knoedlseder + */ + +#ifndef GRMF_HPP +#define GRMF_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GBase.hpp" +#include "GEbounds.hpp" +#include "GMatrixSparse.hpp" +#include "GFits.hpp" +#include "GFitsTable.hpp" + + +/***********************************************************************//** + * @class GRmf + * + * @brief Redistribution Matrix File class + * + * Matrix rows are true energy, columns are channels + ***************************************************************************/ +class GRmf : public GBase { + +public: + // Constructors and destructors + GRmf(void); + explicit GRmf(const std::string& filename); + explicit GRmf(const GEbounds& etrue, const GEbounds& emeasured); + GRmf(const GRmf& rmf); + virtual ~GRmf(void); + + // Operators + GRmf& operator=(const GRmf& rmf); + double& operator()(const int& itrue, const int& imeasured); + const double& operator()(const int& itrue, const int& imeasured) const; + + // Methods + void clear(void); + GRmf* clone(void) const; + int size(void) const; + int ntrue(void) const; + int nmeasured(void) const; + double& at(const int& itrue, const int& imeasured); + const double& at(const int& itrue, const int& imeasured) const; + 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 write(GFits& fits) const; + const std::string& filename(void) const; + std::string print(const GChatter& chatter = NORMAL) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GRmf& rmf); + void free_members(void); + + // Protected members + mutable std::string m_filename; //!< Filename of origin + GEbounds m_ebds_true; //!< True energy boundaries + GEbounds m_ebds_measured; //!< Measured energy boundaries + GMatrixSparse m_matrix; //!< Sparse redistribution matrix + +}; + + +/***********************************************************************//** + * @brief Return content of redistribution matrix bin + * + * @param[in] itrue True energy index [0,...,ntrue()-1]. + * @param[in] imeasured Measured energy index [0,...,nmeasured()-1]. + * + * Returns reference to content of redistribution matrix bin bin with true + * energy index @p itrue and measured energy index @p imeasured. + ***************************************************************************/ +inline +double& GRmf::operator()(const int& itrue, const int& imeasured) +{ + return (m_matrix(itrue, imeasured)); +} + + +/***********************************************************************//** + * @brief Return content of redistribution matrix bin (const version) + * + * @param[in] itrue True energy index [0,...,ntrue()-1]. + * @param[in] imeasured Measured energy index [0,...,nmeasured()-1]. + * + * Returns reference to content of redistribution matrix bin bin with true + * energy index @p itrue and measured energy index @p imeasured. + ***************************************************************************/ +inline +const double& GRmf::operator()(const int& itrue, const int& imeasured) const +{ + return (m_matrix(itrue, imeasured)); +} + + +/***********************************************************************//** + * @brief Return number of redistribution matrix bins + * + * @return Number of redistribution matrix bins. + * + * Returns the number of redistribution matrix bins. + ***************************************************************************/ +inline +int GRmf::size(void) const +{ + return (m_matrix.rows()*m_matrix.columns()); +} + + +/***********************************************************************//** + * @brief Return number of true energy bins in redistribution matrix + * + * @return Number true energy bins in redistribution matrix. + * + * Returns the number of true energy bins in redistribution matrix. + ***************************************************************************/ +inline +int GRmf::ntrue(void) const +{ + return (m_matrix.rows()); +} + + +/***********************************************************************//** + * @brief Return number of measured energy bins in redistribution matrix + * + * @return Number measured energy bins in redistribution matrix. + * + * Returns the number of measured energy bins in redistribution matrix. + ***************************************************************************/ +inline +int GRmf::nmeasured(void) const +{ + return (m_matrix.columns()); +} + + +/***********************************************************************//** + * @brief Return true energy boundaries + * + * @return True energy boundaries for redistribution matrix. + * + * Returns the true energy boundaries for redistribution matrix. + ***************************************************************************/ +inline +const GEbounds& GRmf::etrue(void) const +{ + return m_ebds_true; +} + + +/***********************************************************************//** + * @brief Return measured energy boundaries + * + * @return Measured energy boundaries for redistribution matrix. + * + * Returns the measured energy boundaries for redistribution matrix. + ***************************************************************************/ +inline +const GEbounds& GRmf::emeasured(void) const +{ + return m_ebds_measured; +} + + +/***********************************************************************//** + * @brief Return file name + * + * @return File name from which the RMF information has been read or into + * which RMF information has been saved. + * + * Returns the file name from which the RMF information has been read or into + * which RMF information has been saved. The returned string will be empty if + * no load() or save() method has been called before. + ***************************************************************************/ +inline +const std::string& GRmf::filename(void) const +{ + return (m_filename); +} + +#endif /* GRMF_HPP */ diff --git a/include/GSkyRegion.hpp b/include/GSkyRegion.hpp new file mode 100644 index 0000000..4578c57 --- /dev/null +++ b/include/GSkyRegion.hpp @@ -0,0 +1,134 @@ +/*************************************************************************** + * GSkyRegion.hpp - Abstract virtual sky region base class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GSkyRegion.hpp + * @brief Abstract sky region base class interface definition + * @author Pierrick Martin + */ + +#ifndef GSKYREGION_HPP +#define GSKYREGION_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GBase.hpp" + +/* __ Forward declarations _______________________________________________ */ +class GSkyPixel; +class GSkyDir; + +/***********************************************************************//** + * @class GSkyRegion + * + * @brief Abstract interface for the sky region class + * + * This class provides an abstract interface for a sky region. The sky region + * is defined by an array of parameters names and values specific to the + * derived class where the region type or shape is implemented. + * + * Accessible information elements are: + * - Type of the region (circle, rectangle,...) + * - Solid angle subtended by the region + * + * The object can be initialised from a string in the DS9 region file format, + * and return its description as a string in the DS9 region file format. The + * input/ouput to/from files is handled in the container class GSkyRegions. + ***************************************************************************/ +class GSkyRegion : public GBase { + +public: + // Constructors and destructors + GSkyRegion(void); + GSkyRegion(const GSkyRegion& region); + virtual ~GSkyRegion(void); + + // Operators + virtual GSkyRegion& operator=(const GSkyRegion& region); + + // Pure virtual methods + virtual void clear(void) = 0; + virtual GSkyRegion* clone(void) const = 0; + virtual void read(const std::string& regstring) = 0; + virtual std::string write(void) const = 0; + virtual bool contains(const GSkyDir& dir) const = 0; + virtual bool overlaps(const GSkyRegion& reg) const = 0; + virtual bool contains(const GSkyRegion& reg) const = 0; + virtual std::string print(const GChatter& chatter = NORMAL) const = 0; + + // Implemented methods + std::string type(void) const; + std::string name(void) const; + double solidangle(void) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GSkyRegion& region); + void free_members(void); + + // Protected members + std::string m_type; //!< Type of the region (circle, rectangle,...) + std::string m_name; //!< Name of the region + double m_solid; //!< Solid angle subtended by the region +}; + + +/***********************************************************************//** + * @brief Return region name + * + * @return region name + * + * Returns the region name. + ***************************************************************************/ +inline +std::string GSkyRegion::name(void) const +{ + return (m_name); +} + +/***********************************************************************//** + * @brief Return region type + * + * @return region type + * + * Returns the region type. + ***************************************************************************/ +inline +std::string GSkyRegion::type(void) const +{ + return (m_type); +} + +/***********************************************************************//** + * @brief Return solid angle + * + * @return Solid angle + * + * Returns the solid angle subtended by the region (in sr). + ***************************************************************************/ +inline +double GSkyRegion::solidangle(void) const +{ + return (m_solid); +} + +#endif /* GSKYREGION_HPP */ diff --git a/include/GSkyRegionCircle.hpp b/include/GSkyRegionCircle.hpp new file mode 100644 index 0000000..96b6784 --- /dev/null +++ b/include/GSkyRegionCircle.hpp @@ -0,0 +1,161 @@ +/*************************************************************************** + * GSkyRegionCircle.hpp - circular sky region class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GSkyRegionCircle.hpp + * @brief Abstract sky region base class interface definition + * @author Michael Mayer + */ + +#ifndef GSKYREGIONCIRCLE_HPP +#define GSKYREGIONCIRCLE_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GMath.hpp" +#include "GSkyDir.hpp" +#include "GSkyRegion.hpp" + + +/***********************************************************************//** + * @class GSkyRegionCircle + * + * @brief Abstract interface for the circular sky region class + * + * This class provides an implementation for a circular sky region. The sky + * region is defined by an array of parameters, the meaning of which is + * specific to the derived class where the region type or shape is defined. + * The parameters are GModelPar objects for convenience. + * + * The class holds several properties such as solid angle subtended by the + * region and computed through internal method compute_solid(). + * + * To be clarified: + * - Do we want a member relating the region to an observation run ? + * - Constructor and read/write using XML may not be needed if we use DS9 + * region file format ? + * - Replace GModelPar by double for the parameters (GModelPar is overkill) ? + * + ***************************************************************************/ +class GSkyRegionCircle : public GSkyRegion { + +public: + // Constructors and destructors + GSkyRegionCircle(void); + GSkyRegionCircle(const GSkyRegionCircle& circle); + explicit GSkyRegionCircle(const std::string &line); + GSkyRegionCircle(GSkyDir ¢re, const double &radius); + // GSkyRegionCircle(const double &ra, const double &dec, const double &radius); + virtual ~GSkyRegionCircle(void); + + // Operators + GSkyRegionCircle& operator=(const GSkyRegionCircle& region); + + // Implemented methods + void clear(void); + GSkyRegionCircle* clone(void) const; + GSkyDir centre(void) const; + double radius(void) const; + void radius(const double& radius); + void centre(const GSkyDir& centre); + void centre(const double& ra,const double& dec); + void read(const std::string& line); + std::string write(void) const; + bool contains(const GSkyDir& dir) const; + bool contains(const GSkyRegion& reg) const; + bool overlaps(const GSkyRegion& reg) const; + std::string print(const GChatter& chatter = NORMAL) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GSkyRegionCircle& region); + void free_members(void); + void compute_solid(void); + + // Protected members + GSkyDir m_centre; //!< Centre or reference point of the region + double m_radius; //!< Radius of circular the region [deg] +}; + + +/***********************************************************************//** + * @brief Return region type + * + * @return region type + * + * Returns the region type. + ***************************************************************************/ +inline +double GSkyRegionCircle::radius(void) const +{ + return (m_radius); +} + +/***********************************************************************//** + * @brief Return region centre + * + * @return region centre + * + * Returns the region centre. + ***************************************************************************/ +inline +GSkyDir GSkyRegionCircle::centre(void) const +{ + return (m_centre); +} + + + +/***********************************************************************//** + * @brief Set centre of region + * + * @param[in] dir Center of region. + ***************************************************************************/ +inline +void GSkyRegionCircle::centre(const GSkyDir& dir) +{ + // Set centre + m_centre = dir; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set centre values + * + * @param[in] ra Right ascension value. + * @param[in] dec declination value. + ***************************************************************************/ +inline +void GSkyRegionCircle::centre(const double& ra, const double& dec) +{ + // Set centre values + m_centre.radec_deg(ra,dec); + + // Return + return; +} + +#endif /* GSKYREGIONCIRCLE_HPP */ + diff --git a/include/GSkyRegions.hpp b/include/GSkyRegions.hpp new file mode 100644 index 0000000..62b1bf1 --- /dev/null +++ b/include/GSkyRegions.hpp @@ -0,0 +1,193 @@ +/*************************************************************************** + * GSkyRegions.hpp - Sky region container class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GSkyRegions.hpp + * @brief Sky regions container class definition + * @author Pierrick Martin + */ + +#ifndef GSKYREGIONS_HPP +#define GSKYREGIONS_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GBase.hpp" +#include "GSkyRegion.hpp" +#include "GException.hpp" +#include "GSkyDir.hpp" + +/***********************************************************************//** + * @class GSkyRegions + * + * @brief Sky region container class + * + * This container class collects sky regions objects derived from the + * GSkyRegion abstract class. These can be accessed from the access + * operator[] or through the at() method with index checking. Access to + * the region can be done from index or name, so the name of each region + * in the container has to be unique. + * + * The object can be initialised from a DS9 region file, and has load and + * save methods from/to a DS9 region file. + * + ***************************************************************************/ +class GSkyRegions : public GBase { + +public: + // Constructors and destructors + GSkyRegions(void); + GSkyRegions(const GSkyRegions& regions); + explicit GSkyRegions(const std::string& filename); + virtual ~GSkyRegions(void); + + // Operators + GSkyRegions& operator=(const GSkyRegions& regions); + GSkyRegion* operator[](const int& index); + const GSkyRegion* operator[](const int& index) const; + GSkyRegion* operator[](const std::string& name); + const GSkyRegion* operator[](const std::string& name) const; + + // Methods + void clear(void); + GSkyRegions* clone(void) const; + GSkyRegion* at(const int& index); + const GSkyRegion* at(const int& index) const; + int size(void) const; + bool isempty(void) const; + GSkyRegion* set(const int& index, const GSkyRegion& region); + GSkyRegion* set(const std::string& name, const GSkyRegion& region); + GSkyRegion* append(const GSkyRegion& region); + GSkyRegion* insert(const int& index, const GSkyRegion& region); + GSkyRegion* insert(const std::string& name, const GSkyRegion& region); + void remove(const int& index); + void remove(const std::string& name); + void reserve(const int& num); + void extend(const GSkyRegions& regions); + bool contains(const std::string& name) const; + bool contains(const GSkyDir& dir) const; + bool overlaps(const GSkyRegion& reg) const; + void load(const std::string& filename); + void save(const std::string& filename) const; + const std::string& filename(void) const; + std::string print(const GChatter& chatter = NORMAL) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GSkyRegions& regions); + void free_members(void); + int get_index(const std::string& name) const; + + // Protected members + mutable std::string m_filename; //!< Filename of origin + std::vector m_regions; //!< List of regions +}; + + +/***********************************************************************//** + * @brief Return pointer to region + * + * @param[in] index region index [0,...,size()-1]. + * + * Returns a pointer to the region with the specified @p index. + ***************************************************************************/ +inline +GSkyRegion* GSkyRegions::operator[](const int& index) +{ + return (m_regions[index]); +} + + +/***********************************************************************//** + * @brief Return pointer to region (const version) + * + * @param[in] index region index [0,...,size()-1]. + * + * Returns a const pointer to the region with the specified @p index. + ***************************************************************************/ +inline +const GSkyRegion* GSkyRegions::operator[](const int& index) const +{ + return (m_regions[index]); +} + + +/***********************************************************************//** + * @brief Return number of regions in container + * + * @return Number of regions in container. + * + * Returns the number of regions in the region container. + ***************************************************************************/ +inline +int GSkyRegions::size(void) const +{ + return (m_regions.size()); +} + + +/***********************************************************************//** + * @brief Signals if there are no regions in container + * + * @return True if container is empty, false otherwise. + * + * Signals if the region container does not contain any region. + ***************************************************************************/ +inline +bool GSkyRegions::isempty(void) const +{ + return (m_regions.empty()); +} + + +/***********************************************************************//** + * @brief Reserves space for regions in container + * + * @param[in] num Number of regions + * + * Reserves space for @p num regions in the container. + ***************************************************************************/ +inline +void GSkyRegions::reserve(const int& num) +{ + m_regions.reserve(num); + return; +} + + +/***********************************************************************//** + * @brief Return file name + * + * @return File name from which the regions have been loaded or into which + * the regions have been saved. + * + * Returns the file name from which the regions have been loaded or into + * which the regions have been saved. The returned string will be empty if + * no load() or save() method has been called before. + ***************************************************************************/ +inline +const std::string& GSkyRegions::filename(void) const +{ + return (m_filename); +} + +#endif /* GSKYREGIONS_HPP */ diff --git a/include/GTools.hpp b/include/GTools.hpp index 1e05e71..48026cf 100644 --- a/include/GTools.hpp +++ b/include/GTools.hpp @@ -87,6 +87,7 @@ namespace gammalib { bool file_exists(const std::string& filename); bool isinfinite(const double& x); bool isnotanumber(const double& x); + bool contains(std::string str1, std::string str2); } /***********************************************************************//** diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 6427079..b5e61da 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -107,6 +107,11 @@ /* __ VO module __________________________________________________________ */ #include "GVOClient.hpp" +/* __ XSPEC module _______________________________________________________ */ +#include "GArf.hpp" +#include "GPha.hpp" +#include "GRmf.hpp" + /* __ Application module _________________________________________________ */ #include "GApplication.hpp" #include "GLog.hpp" @@ -132,6 +137,9 @@ #include "GSkyDir.hpp" #include "GSkyPixel.hpp" #include "GSkymap.hpp" +#include "GSkyRegions.hpp" +#include "GSkyRegion.hpp" +#include "GSkyRegionCircle.hpp" #include "GWcs.hpp" #include "GWcsRegistry.hpp" #include "GWcslib.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index fdf1d66..01a6eaa 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -61,6 +61,9 @@ pkginclude_HEADERS = GammaLib.hpp \ GXmlAttribute.hpp \ GXmlPI.hpp \ GVOClient.hpp \ + GPha.hpp \ + GArf.hpp \ + GRmf.hpp \ GApplication.hpp \ GLog.hpp \ GPars.hpp \ @@ -75,6 +78,9 @@ pkginclude_HEADERS = GammaLib.hpp \ GSkyDir.hpp \ GSkyPixel.hpp \ GSkymap.hpp \ + GSkyRegion.hpp \ + GSkyRegions.hpp \ + GSkyRegionCircle.hpp \ GWcs.hpp \ GWcsRegistry.hpp \ GWcslib.hpp \ diff --git a/inst/cta/Makefile.am b/inst/cta/Makefile.am index d4b7524..6776831 100644 --- a/inst/cta/Makefile.am +++ b/inst/cta/Makefile.am @@ -11,6 +11,8 @@ INCLUDES = -I$(top_srcdir)/include \ sources = src/GCTASupport.cpp \ src/GCTAException.cpp \ src/GCTAObservation.cpp \ + src/GCTAOnOffObservation.cpp \ + src/GCTAOnOffObservations.cpp \ src/GCTAEventList.cpp \ src/GCTAEventAtom.cpp \ src/GCTAEventCube.cpp \ @@ -41,6 +43,8 @@ sources = src/GCTASupport.cpp \ # Define headers to be installed pkginclude_HEADERS = include/GCTAException.hpp \ include/GCTAObservation.hpp \ + include/GCTAOnOffObservation.hpp \ + include/GCTAOnOffObservations.hpp \ include/GCTAEventList.hpp \ include/GCTAEventAtom.hpp \ include/GCTAEventCube.hpp \ diff --git a/inst/cta/include/GCTAOnOffObservation.hpp b/inst/cta/include/GCTAOnOffObservation.hpp new file mode 100644 index 0000000..0b58f5f --- /dev/null +++ b/inst/cta/include/GCTAOnOffObservation.hpp @@ -0,0 +1,102 @@ +/*************************************************************************** + * GCTAOnOffObservation.hpp - CTA on-off observation class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Chia-Chun Lu & Christoph Deil * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAOnOffObservation.hpp + * @brief CTA on-off observation class definition + * @author Chia-Chun Lu & Christoph Deil + */ + +#ifndef GCTAONOFFOBSERVATION_HPP +#define GCTAONOFFOBSERVATION_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GBase.hpp" +#include "GPha.hpp" +#include "GArf.hpp" +#include "GRmf.hpp" +#include "GCTAEventList.hpp" +#include "GCTAEventAtom.hpp" +#include "GCTAObservation.hpp" +#include "GSkyRegions.hpp" + + +/***********************************************************************//** + * @class GCTAOnOffObservation + * + * @brief CTA on-off observation class + ***************************************************************************/ +class GCTAOnOffObservation : public GBase { + +public: + // Constructors and destructors + GCTAOnOffObservation(void); + GCTAOnOffObservation(const GEbounds& ereco, const GSkyRegions& on, + const GSkyRegions& off); + GCTAOnOffObservation(const GCTAOnOffObservation& obs); + virtual ~GCTAOnOffObservation(void); + + // Operators + GCTAOnOffObservation& operator=(const GCTAOnOffObservation& obs); + + // Methods + void clear(void); + GCTAOnOffObservation* clone(void) const; + void name(const std::string& name) { m_name = name; } + void instrument(const std::string& instrument) { m_instrument = instrument; } + void id(const std::string& id) { m_id = id; } + void off_regions(const GSkyRegions& regions) {m_off_regions = regions;} + void on_regions(const GSkyRegions& regions) {m_on_regions = regions;} + const std::string& name(void) const { return m_name; } + const std::string& instrument(void) const { return m_instrument; } + const std::string& id(void) const { return m_id; } + const GPha& on_spec(void) const { return m_on_spec; } + const GPha& off_spec(void) const { return m_off_spec; } + const GArf& arf(void) { return m_arf; }; + const GRmf& rmf(void) {return m_rmf;} + void fill(const GCTAObservation& obs); + void compute_response(const GCTAObservation& obs, + const GEbounds& etrue); + void read(const GXmlElement& xml); + void write(GXmlElement& xml) const; + std::string print(const GChatter& chatter = NORMAL) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GCTAOnOffObservation& obs); + void free_members(void); + void compute_arf(const GCTAObservation& obs); + void compute_rmf(const GCTAObservation& obs, const GEbounds& etrue); + + // Protected data members + std::string m_name; //!< Name + std::string m_instrument; //!< Instrument name + std::string m_id; //!< Observation identifier + GPha m_on_spec; + GPha m_off_spec; + GArf m_arf; + GRmf m_rmf; + GSkyRegions m_on_regions; + GSkyRegions m_off_regions; +}; + +#endif /* GCTAONOFFOBSERVATION_HPP */ diff --git a/inst/cta/include/GCTAOnOffObservations.hpp b/inst/cta/include/GCTAOnOffObservations.hpp new file mode 100644 index 0000000..57ef8ac --- /dev/null +++ b/inst/cta/include/GCTAOnOffObservations.hpp @@ -0,0 +1,209 @@ +/*************************************************************************** + * GCTAOnOffObservations.hpp - CTA on-off observation container class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAOnOffObservations.hpp + * @brief CTA on-off observation container class definition + * @author Pierrick Martin + */ + +#ifndef GCTAONOFFOBSERVATIONS_HPP +#define GCTAONOFFOBSERVATIONS_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include +#include "GContainer.hpp" +#include "GCTAOnOffObservation.hpp" +#include "GModels.hpp" + + +/***********************************************************************//** + * @class GCTAOnOffObservations + * + * @brief ON/OFF Observation container class + * + * This class in a container of GCTAOnOffObservation objects. Still some + * work to do and things to be clarified... + ***************************************************************************/ +class GCTAOnOffObservations : public GContainer { + +public: + // Constructors and destructors + GCTAOnOffObservations(void); + GCTAOnOffObservations(const GCTAOnOffObservations& obs); + explicit GCTAOnOffObservations(const std::string& filename); + virtual ~GCTAOnOffObservations(void); + + // Operators + GCTAOnOffObservations& operator=(const GCTAOnOffObservations& obs); + GCTAOnOffObservation* operator[](const int& index); + const GCTAOnOffObservation* operator[](const int& index) const; + + // Methods + void clear(void); + GCTAOnOffObservations* clone(void) const; + GCTAOnOffObservation* at(const int& index); + const GCTAOnOffObservation* at(const int& index) const; + int size(void) const; + bool isempty(void) const; + GCTAOnOffObservation* set(const int& index, const GCTAOnOffObservation& obs); + GCTAOnOffObservation* append(const GCTAOnOffObservation& obs); + GCTAOnOffObservation* insert(const int& index, const GCTAOnOffObservation& obs); + void remove(const int& index); + void reserve(const int& num); + void extend(const GCTAOnOffObservations& obs); + bool contains(const std::string& instrument, + const std::string& id) const; + void load(const std::string& filename); + void save(const std::string& filename) const; + void read(const GXml& xml); + void write(GXml& xml) const; + void models(const GModels& models); + void models(const std::string& filename); + const GModels& models(void) const; + std::string print(const GChatter& chatter = NORMAL) const; + + /* + // To do + Implement fit functions for both cases summed over runs or not ? + Implement significance or have it at ctool level ? + Implement flux points computation and fit of a spectral function ? + // Do we want that ? If yes, need to retrieve inline functions from GObservations... + double npred(void) const; + // Left aside at the moment + void optimize(GOptimizer& opt); + */ + +protected: + // Protected methods + void init_members(void); + void copy_members(const GCTAOnOffObservations& obs); + void free_members(void); + int get_index(const std::string& instrument, + const std::string& id) const; + + // Protected members + std::vector m_obs; //!< List of observations + GModels m_models; //!< List of models +}; + + +/***********************************************************************//** + * @brief Return pointer to observation + * + * @param[in] index Observation index [0,...,size()-1]. + * @return Observation. + * + * Returns a pointer to the observation with specified @p index. + ***************************************************************************/ +inline +GCTAOnOffObservation* GCTAOnOffObservations::operator[](const int& index) +{ + return m_obs[index]; +} + + +/***********************************************************************//** + * @brief Return pointer to observation (const version) + * + * @param[in] index Observation index [0,...,size()-1]. + * + * Returns a const pointer to the observation with specified @p index. + ***************************************************************************/ +inline +const GCTAOnOffObservation* GCTAOnOffObservations::operator[](const int& index) const +{ + return m_obs[index]; +} + + +/***********************************************************************//** + * @brief Return number of observations in container + * + * @return Number of observations in container. + * + * Returns the number of observations in the observation container. + ***************************************************************************/ +inline +int GCTAOnOffObservations::size(void) const +{ + return (m_obs.size()); +} + + +/***********************************************************************//** + * @brief Signals if there are no observations in container + * + * @return True if container is empty, false otherwise. + * + * Signals if the observation container does not contain any observation. + ***************************************************************************/ +inline +bool GCTAOnOffObservations::isempty(void) const +{ + return (m_obs.empty()); +} + + +/***********************************************************************//** + * @brief Reserves space for observations in container + * + * @param[in] num Number of observations. + * + * Reserves space for @p num observations in the container. + ***************************************************************************/ +inline +void GCTAOnOffObservations::reserve(const int& num) +{ + m_obs.reserve(num); + return; +} + + +/***********************************************************************//** + * @brief Set model container + * + * @param[in] models Model container. + * + * Sets the model container for the observations. + ***************************************************************************/ +inline +void GCTAOnOffObservations::models(const GModels& models) +{ + m_models = models; + return; +} + + +/***********************************************************************//** + * @brief Return model container + * + * @return Model container. + * + * Returns the model container of the observations. + ***************************************************************************/ +inline +const GModels& GCTAOnOffObservations::models(void) const +{ + return m_models; +} + +#endif /* GCTAONOFFOBSERVATIONS_HPP */ diff --git a/inst/cta/pyext/GCTAOnOffObservation.i b/inst/cta/pyext/GCTAOnOffObservation.i new file mode 100644 index 0000000..1b92d3d --- /dev/null +++ b/inst/cta/pyext/GCTAOnOffObservation.i @@ -0,0 +1,76 @@ +/*************************************************************************** + * GCTAOnOffObservation.i - CTA on-off observation class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAOnOffObservation.i + * @brief CTA on-off observation class definition + * @author Michael Mayer + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GCTAOnOffObservation.hpp" +%} + + +/***********************************************************************//** + * @class GCTAOnOffObservation + * + * @brief CTA on-off observation class + ***************************************************************************/ +class GCTAOnOffObservation : public GBase { +public: + // Constructors and destructors + GCTAOnOffObservation(void); + GCTAOnOffObservation(const GEbounds& ereco, const GSkyRegions& on, + const GSkyRegions& off); + GCTAOnOffObservation(const GCTAOnOffObservation& obs); + virtual ~GCTAOnOffObservation(void); + + // Methods + void clear(void); + GCTAOnOffObservation* clone(void) const; + void name(const std::string& name) { m_name = name; } + void instrument(const std::string& instrument) { m_instrument = instrument; } + void id(const std::string& id) { m_id = id; } + void off_regions(const GSkyRegions& regions) {m_off_regions = regions;} + void on_regions(const GSkyRegions& regions) {m_on_regions = regions;} + const std::string& name(void) const { return m_name; } + const std::string& instrument(void) const { return m_instrument; } + const std::string& id(void) const { return m_id; } + const GPha& on_spec(void) const { return m_on_spec; } + const GPha& off_spec(void) const { return m_off_spec; } + const GArf& arf(void) { return m_arf; }; + const GRmf& rmf(void) {return m_rmf;} + void fill(const GCTAObservation& obs); + void compute_response(const GCTAObservation& obs, + const GEbounds& etrue); + void read(const GXmlElement& xml); + void write(GXmlElement& xml) const; +}; + + +/***********************************************************************//** + * @brief GCTOnOffAObservation class extension + ***************************************************************************/ +%extend GCTAOnOffObservation { + GCTAOnOffObservation copy() { + return (*self); + } +} diff --git a/inst/cta/pyext/GCTAOnOffObservations.i b/inst/cta/pyext/GCTAOnOffObservations.i new file mode 100644 index 0000000..fc7cc8d --- /dev/null +++ b/inst/cta/pyext/GCTAOnOffObservations.i @@ -0,0 +1,97 @@ +/*************************************************************************** + * GCTAOnOffObservations.i - CTA on-off observation container class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAOnOffObservations.i + * @brief CTA on-off observation container class definition + * @author Pierrick Martin + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GCTAOnOffObservations.hpp" +%} + + +/***********************************************************************//** + * @class GCTAOnOffObservations + * + * @brief ON/OFF Observation container class + ***************************************************************************/ +class GCTAOnOffObservations : public GContainer { +public: + // Constructors and destructors + GCTAOnOffObservations(void); + GCTAOnOffObservations(const GCTAOnOffObservations& obs); + explicit GCTAOnOffObservations(const std::string& filename); + virtual ~GCTAOnOffObservations(void); + + // Methods + void clear(void); + GCTAOnOffObservations* clone(void) const; + GCTAOnOffObservation* at(const int& index); + const GCTAOnOffObservation* at(const int& index) const; + int size(void) const; + bool isempty(void) const; + GCTAOnOffObservation* set(const int& index, const GCTAOnOffObservation& obs); + GCTAOnOffObservation* append(const GCTAOnOffObservation& obs); + GCTAOnOffObservation* insert(const int& index, const GCTAOnOffObservation& obs); + void remove(const int& index); + void reserve(const int& num); + void extend(const GCTAOnOffObservations& obs); + bool contains(const std::string& instrument, + const std::string& id) const; + void load(const std::string& filename); + void save(const std::string& filename) const; + void read(const GXml& xml); + void write(GXml& xml) const; + void models(const GModels& models); + void models(const std::string& filename); + const GModels& models(void) const; +}; + + +/***********************************************************************//** + * @brief GCTAOnOffObservations class extension + ***************************************************************************/ +%extend GCTAOnOffObservations { + GCTAOnOffObservation* __getitem__(const int& index) { + if (index >= 0 && index < self->size()) { + return (*self)[index]; + } + else { + throw GException::out_of_range("__getitem__(int)", index, self->size()); + } + } + void __setitem__(const int& index, const GCTAOnOffObservation& val) { + if (index>=0 && index < self->size()) { + self->set(index, val); + return; + } + else { + throw GException::out_of_range("__setitem__(int)", index, self->size()); + } + } + int __len__() { + return (self->size()); + } + GCTAOnOffObservations copy() { + return (*self); + } +}; diff --git a/inst/cta/pyext/cta.i b/inst/cta/pyext/cta.i index 0fdafe5..ce16010 100644 --- a/inst/cta/pyext/cta.i +++ b/inst/cta/pyext/cta.i @@ -63,6 +63,8 @@ /* __ CTA ________________________________________________________________ */ %include "GCTAObservation.i" +%include "GCTAOnOffObservation.i" +%include "GCTAOnOffObservations.i" %include "GCTAEventCube.i" %include "GCTAEventList.i" %include "GCTAEventBin.i" diff --git a/inst/cta/src/GCTAOnOffObservation.cpp b/inst/cta/src/GCTAOnOffObservation.cpp new file mode 100644 index 0000000..cc5ff5c --- /dev/null +++ b/inst/cta/src/GCTAOnOffObservation.cpp @@ -0,0 +1,642 @@ +/*************************************************************************** + * GCTAOnOffObservation.cpp - CTA on-off observation class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Chia-Chun Lu & Christoph Deil * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAOnOffObservation.cpp + * @brief CTA on-off observation class implementation + * @author Chia-Chun Lu & Christoph Deil + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GCTAOnOffObservation.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_WRITE "GCTAOnOffObservation::write(GXmlElement&)" +#define G_READ "GCTAOnOffObservation::read(GXmlElement&)" +#define G_FILL "GCTAOnOffObservation::fill(GCTAObservation&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GCTAOnOffObservation::GCTAOnOffObservation(void) +{ + // Initialise private members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] obs CTA on-off observation. + ***************************************************************************/ +GCTAOnOffObservation::GCTAOnOffObservation(const GCTAOnOffObservation& obs) +{ + // Initialise private + init_members(); + + // Copy members + copy_members(obs); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Parameter constructor + * + * @param[in] ereco Reconstructed energy bins. + * @param[in] on On regions. + * @param[in] off Off regions. + ***************************************************************************/ +GCTAOnOffObservation::GCTAOnOffObservation(const GEbounds& ereco, + const GSkyRegions& on, + const GSkyRegions& off) +{ + // Initialise private + init_members(); + + // Initialise spectra + m_on_spec = GPha(ereco); + m_off_spec = GPha(ereco); + + // Store regions + m_on_regions = on; + m_off_regions = off; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GCTAOnOffObservation::~GCTAOnOffObservation(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] obs CTA on-off observation. + ***************************************************************************/ +GCTAOnOffObservation& GCTAOnOffObservation::operator=(const GCTAOnOffObservation& obs) +{ + // Execute only if object is not identical + if (this != &obs) { + + // Free members + free_members(); + + // Initialise private members for clean destruction + init_members(); + + // Copy members + copy_members(obs); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear instance + * + * This method properly resets the object to an initial state. + ***************************************************************************/ +void GCTAOnOffObservation::clear(void) +{ + // Free class members + free_members(); + + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone instance + **************************************************************************/ +GCTAOnOffObservation* GCTAOnOffObservation::clone(void) const +{ + return new GCTAOnOffObservation(*this); +} + + +/***********************************************************************//** + * @brief read observation from an xml element + * + * @param[in] xml GXmlElement. + ***************************************************************************/ +void GCTAOnOffObservation::read(const GXmlElement& xml) +{ + // clean object + clear(); + + // Extract instrument name + m_instrument = xml.attribute("instrument"); + + // Determine number of parameter nodes in XML element + int npars = xml.elements("parameter"); + + // Verify that XML element has exactly 6 parameters + if (xml.elements() != 6 || npars != 6) { + throw GException::xml_invalid_parnum(G_READ, xml, + "CTA observation requires exactly 6 parameters."); + } + + // Extract parameters + int npar[] = {0, 0, 0, 0, 0, 0}; + for (int i = 0; i < npars; ++i) { + + // Get parameter element + const GXmlElement* par = xml.element("parameter", i); + + // Handle Pha_on + if (par->attribute("name") == "Pha_on") { + + // Read pha_on file name + std::string filename = par->attribute("file"); + + // Load pha_on + m_on_spec.load(filename); + + // Increment parameter counter + npar[0]++; + } + + // Handle Pha_off + else if (par->attribute("name") == "Pha_off") { + + // Read countsmap file name + std::string filename = par->attribute("file"); + + // Load pha_off + m_off_spec.load(filename); + + // Increment parameter counter + npar[1]++; + } + + // Handle on regions + else if(par->attribute("name") == "Regions_on") { + + // Get filename + std::string filename = par->attribute("file"); + + // load on regions + m_on_regions.load(filename); + + // Increase number of parameters + npar[2]++; + } + + // Handle off regions + else if(par->attribute("name") == "Regions_off") { + + // Get filename + std::string filename = par->attribute("file"); + + // load off regions + m_on_regions.load(filename); + + // Increase number of parameters + npar[3]++; + } + + // Handle Arf + else if(par->attribute("name") == "Arf") { + + // Get filename + std::string filename = par->attribute("file"); + + // load arf + m_arf.load(filename); + + // Increase number of parameters + npar[4]++; + } + + // Handle Rmf + else if(par->attribute("name") == "Rmf") { + + // Get filename + std::string filename = par->attribute("file"); + + // load Rmf + m_rmf.load(filename); + + // Increase number of parameters + npar[5]++; + } + + } // endfor: looped over all parameters + + // Verify that all parameters were found + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1 || npar[4] != 1 || npar[5] != 1) { + throw GException::xml_invalid_parnames(G_READ, xml, + "Require \"Pha_on\" or \"Pha_off\" and \"Regions_on\"" + ", \"Regions_off\",\"Rmf\" and \"Arf\" parameters."); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief write observation to an xml element + * + * @param[in] xml GXmlElement. + ***************************************************************************/ +void GCTAOnOffObservation::write(GXmlElement& xml) const +{ + // If XML element has 0 nodes then append 6 parameter nodes + if (xml.elements() == 0) { + xml.append(GXmlElement("parameter name=\"Pha_on\"")); + xml.append(GXmlElement("parameter name=\"Pha_off\"")); + xml.append(GXmlElement("parameter name=\"Regions_on\"")); + xml.append(GXmlElement("parameter name=\"Regions_off\"")); + xml.append(GXmlElement("parameter name=\"Arf\"")); + xml.append(GXmlElement("parameter name=\"Rmf\"")); + } + + // Verify that XML element has exactly 4 parameters + if (xml.elements() != 6 || xml.elements("parameter") != 6) { + throw GException::xml_invalid_parnum(G_WRITE, xml, + "CTAOnOffObservation requires exactly 6 parameters."); + } + + // Set or update parameter attributes + int npar[] = {0, 0, 0, 0, 0, 0}; + for (int i = 0; i < 6; ++i) { + + // Get parameter element + GXmlElement* par = xml.element("parameter", i); + + // Handle on counts + if (par->attribute("name") == "Pha_on") { + par->attribute("file", m_on_spec.filename()); + npar[0]++; + } + + // Handle off counts + else if (par->attribute("name") == "Pha_off") { + par->attribute("file", m_off_spec.filename()); + npar[0]++; + } + + // Handle on regions + else if (par->attribute("name") == "Regions_on") { + par->attribute("file", m_on_regions.filename()); + npar[0]++; + } + + // Handle off regions + else if (par->attribute("name") == "Regions_off") { + par->attribute("file", m_off_regions.filename()); + npar[0]++; + } + + // Handle effective area + else if (par->attribute("name") == "Arf") { + par->attribute("file", m_arf.filename()); + npar[0]++; + } + + // Handle energy resolution + else if (par->attribute("name") == "Rmf") { + par->attribute("file", m_rmf.filename()); + npar[0]++; + } + } + // Verify that all required parameters are present + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1 || npar[4] != 1 || npar[5] !=1) { + throw GException::xml_invalid_parnames(G_WRITE, xml, + "Require \"Pha_on\" or \"Pha_off\" and \"Regions_on\", \"Regions_off\"" + " \"Arf\", \"Rmf\" parameters."); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Fill events in the on and off spectra + * + * @param[in] obs CTA observation. + * + * @exception GException::invalid_value + * No CTA event list found in CTA observation. + ***************************************************************************/ +void GCTAOnOffObservation::fill(const GCTAObservation& obs) +{ + // Get CTA event list pointer + const GCTAEventList* events = dynamic_cast(obs.events()); + if (events == NULL) { + throw GException::invalid_value(G_FILL, + "GCTAObservation must provide an event list"); + } + + // Loop over all events + for (int i = 0; i < events->size(); ++i) { + + // Get measured event direction + const GCTAEventAtom* atom = (*events)[i]; + GSkyDir dir = atom->dir().dir(); + + // Fill in spectrum according to region containment + if (m_on_regions.contains(dir)) { + m_on_spec.fill(atom->energy()); + } + else if (m_off_regions.contains(dir)) { + m_off_spec.fill(atom->energy()); + } + + } // endfor: looped over all events + + // Return + return; +} + + +/***********************************************************************//** + * @brief Compute response of an on/off observation + * + * @param[in] obs GCTAObservation. + * @param[in] etrue True energy boundaries. + ***************************************************************************/ +void GCTAOnOffObservation::compute_response(const GCTAObservation& obs, + const GEbounds& etrue) +{ + // Compute response components + compute_arf(obs); + compute_rmf(obs, etrue); + + // Return + return; +} + +/***********************************************************************//** + * @brief Compute arf of an on/off observation + * + * @param[in] obs GCTAObservation. + * + * @todo Implement GCTAResponse::npred usage. + ***************************************************************************/ +void GCTAOnOffObservation::compute_arf(const GCTAObservation& obs) +{ + // Set constant response parameters + const double theta = 0.0; + const double phi = 0.0; + const double zenith = 0.0; + const double azimuth = 0.0; + + // Get energy boundaries from on spectrum + GEbounds ereco = m_on_spec.ebounds(); + + // Continue only if there are spectral bins + int nreco = ereco.size(); + if (nreco > 0) { + + // Get CTA response pointer + GCTAResponse* response = obs.response(); + + // Initialize ARF + m_arf = GArf(ereco); + + // Loop over reconstructed energies + for (int i = 0; i < nreco; ++i) { + + // Get mean energy of bin + double logenergy = ereco.elogmean(i).log10TeV(); + + // Set specresp value + m_arf[i] = response->aeff(theta, + phi, + zenith, + azimuth, + logenergy); + + } // endfor: looped over reconstructed energies + + } // endif: there were energy bins + + // Return + return; +} + + +/***********************************************************************//** + * @brief compute rmf of an on/off observation + * + * @param[in] obs GCTAObservation. + * @param[in] etrue True energy boundaries. + ***************************************************************************/ +void GCTAOnOffObservation::compute_rmf(const GCTAObservation& obs, + const GEbounds& etrue) +{ + // Set constant response parameters + const double theta = 0.0; + const double phi = 0.0; + const double zenith = 0.0; + const double azimuth = 0.0; + + // Get reconstructed energy boundaries from on spectrum + GEbounds ereco = m_on_spec.ebounds(); + + // Continue only if there are spectral bins + int ntrue = etrue.size(); + int nreco = ereco.size(); + if (ntrue > 0 && nreco > 0) { + + // Get CTA response pointer + GCTAResponse* response = obs.response(); + + // Initialize RMF + m_rmf = GRmf(etrue, ereco); + + // Loop over reconstructed energy + for (int ireco = 0; ireco < nreco; ++ireco) { + + // Compute reconstructed energy + double eng_reco = ereco.elogmean(ireco).log10TeV(); + + // Loop over true energy + for (int itrue = 0; itrue < ntrue; ++itrue) { + + // Compute true energy + double eng_true = etrue.elogmean(itrue).log10TeV(); + + // Set RMF value + m_rmf(itrue, ireco) = response->edisp(eng_reco, + theta, + phi, + zenith, + azimuth, + eng_true); + } // endfor: looped over true energy + } // endfor: looped over reconstructed energy + } // endif: there were energy bins + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print class information + ***************************************************************************/ +std::string GCTAOnOffObservation::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GCTAOnOffObservation ==="); + + // Append parameters + result.append("\n"+gammalib::parformat("Name")+m_name); + result.append("\n"+gammalib::parformat("Identifier")+m_id); + + // Append spectra, ARF and RMF + result.append("\n"+m_on_spec.print()); + result.append("\n"+m_off_spec.print()); + result.append("\n"+m_arf.print()); + result.append("\n"+m_rmf.print()); + + // Append regions + result.append("\n"+m_on_regions.print()); + result.append("\n"+m_off_regions.print()); + } + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GCTAOnOffObservation::init_members(void) +{ + // Initialise members + m_name.clear(); + m_id.clear(); + m_instrument = "CTA"; + m_on_spec.clear(); + m_off_spec.clear(); + m_arf.clear(); + m_rmf.clear(); + m_on_regions.clear(); + m_off_regions.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] obs CTA on-off observation. + ***************************************************************************/ +void GCTAOnOffObservation::copy_members(const GCTAOnOffObservation& obs) +{ + // Copy attributes + m_name = obs.m_name; + m_id = obs.m_id; + m_instrument = obs.m_instrument; + m_on_spec = obs.m_on_spec; + m_off_spec = obs.m_off_spec; + m_arf = obs.m_arf; + m_rmf = obs.m_rmf; + m_on_regions = obs.m_on_regions; + m_off_regions = obs.m_off_regions; + + // Return + return; +} + + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GCTAOnOffObservation::free_members(void) +{ + // Return + return; +} diff --git a/inst/cta/src/GCTAOnOffObservations.cpp b/inst/cta/src/GCTAOnOffObservations.cpp new file mode 100644 index 0000000..22023a4 --- /dev/null +++ b/inst/cta/src/GCTAOnOffObservations.cpp @@ -0,0 +1,832 @@ +/*************************************************************************** + * GCTAOnOffObservations.cpp - ON/OFF Observation container class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAOnOffObservations.cpp + * @brief ON/OFF Observations container class implementation + * @author Pierrick Martin + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GTools.hpp" +#include "GException.hpp" +#include "GCTAOnOffObservations.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_AT "GCTAOnOffObservations::at(int&)" +#define G_SET "GCTAOnOffObservations::set(int&, GCTAOnOffObservation&)" +#define G_APPEND "GCTAOnOffObservations::append(GCTAOnOffObservation&)" +#define G_INSERT "GCTAOnOffObservations::insert(int&, GCTAOnOffObservation&)" +#define G_REMOVE "GCTAOnOffObservations::remove(int&)" +#define G_EXTEND "GCTAOnOffObservations::extend(GCTAOnOffObservations&)" +#define G_READ "GCTAOnOffObservations::read(GXml&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GCTAOnOffObservations::GCTAOnOffObservations(void) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param obs Observation container. + ***************************************************************************/ +GCTAOnOffObservations::GCTAOnOffObservations(const GCTAOnOffObservations& obs) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(obs); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Load constructor + * + * @param[in] filename XML filename. + * + * Construct the observation container by loading all relevant information + * from a XML file. Please refer to the read() method for more information + * about the structure of the XML file. + ***************************************************************************/ +GCTAOnOffObservations::GCTAOnOffObservations(const std::string& filename) +{ + // Initialise members + init_members(); + + // Load XML file + load(filename); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GCTAOnOffObservations::~GCTAOnOffObservations(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] obs Observation container. + * @return Observation container. + ***************************************************************************/ +GCTAOnOffObservations& GCTAOnOffObservations::operator=(const GCTAOnOffObservations& obs) +{ + // Execute only if object is not identical + if (this != &obs) { + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(obs); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear observations + ***************************************************************************/ +void GCTAOnOffObservations::clear(void) +{ + // Free members + free_members(); + + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone observations + * + * @return Pointer to deep copy of observation container. + ***************************************************************************/ +GCTAOnOffObservations* GCTAOnOffObservations::clone(void) const +{ + // Clone observations + return new GCTAOnOffObservations(*this); +} + + +/***********************************************************************//** + * @brief Return pointer to observation + * + * @param[in] index Observation index [0,...,size()-1]. + * @return Observation. + * + * @exception GException::out_of_range + * Operation index is out of range. + * + * Returns a pointer to the observation with specified @p index. + ***************************************************************************/ +GCTAOnOffObservation* GCTAOnOffObservations::at(const int& index) +{ + // Raise exception if index is out of range + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_AT, index, 0, size()-1); + } + + // Return pointer + return m_obs[index]; +} + + +/***********************************************************************//** + * @brief Return pointer to observation (const version) + * + * @param[in] index Observation index [0,...,size()-1]. + * + * @exception GException::out_of_range + * Operation index is out of range. + * + * Returns a const pointer to the observation with specified @p index. + ***************************************************************************/ +const GCTAOnOffObservation* GCTAOnOffObservations::at(const int& index) const +{ + // Raise exception if index is out of range + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_AT, index, 0, size()-1); + } + + // Return pointer + return m_obs[index]; +} + + +/***********************************************************************//** + * @brief Set observation in container + * + * @param[in] index Observation index [0,...,size()-1]. + * @param[in] obs Observation. + * @return Pointer to deep copy of observation. + * + * @exception GException::out_of_range + * Observation index is out of range. + * @exception GException::invalid_value + * Observation with same instrument and identifier already + * exists in container. + * + * Set a deep copy and observation @p obs at the specified @p index in the + * container. + ***************************************************************************/ +GCTAOnOffObservation* GCTAOnOffObservations::set(const int& index, const GCTAOnOffObservation& obs) +{ + // Compile option: raise exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_SET, index, 0, size()-1); + } + #endif + + // Raise an exception if an observation with specified instrument and + // identifier already exists + int inx = get_index(obs.instrument(), obs.id()); + if (inx != -1 && inx != index) { + std::string msg = + "Attempt to set \""+obs.instrument()+"\" observation with" + " identifier \""+obs.id()+"\" in observation container at" + " index "+gammalib::str(index)+", but an observation with the" + " same attributes exists already at index "+gammalib::str(inx)+ + " in the container.\n" + "Every observation for a given instrument in the observation" + " container needs a unique identifier."; + throw GException::invalid_value(G_SET, msg); + } + + // Delete existing observation + if (m_obs[index] != NULL) delete m_obs[index]; + + // Store pointer to a deep copy of the observation + m_obs[index] = obs.clone(); + + // Return pointer + return m_obs[index]; +} + + +/***********************************************************************//** + * @brief Append observation to container + * + * @param[in] obs Observation. + * @return Pointer to deep copy of observation. + * + * @exception GException::invalid_value + * Observation with same instrument and identifier already + * exists in container. + * + * Appends a deep copy of an observation to the container. + ***************************************************************************/ +GCTAOnOffObservation* GCTAOnOffObservations::append(const GCTAOnOffObservation& obs) +{ + // Raise an exception if an observation with specified instrument and + // identifier already exists + int inx = get_index(obs.instrument(), obs.id()); + if (inx != -1) { + std::string msg = + "Attempt to append \""+obs.instrument()+"\" observation with" + " identifier \""+obs.id()+"\" to observation container, but an" + " observation with the same attributes exists already at" + " index "+gammalib::str(inx)+" in the container.\n" + "Every observation for a given instrument in the observation" + " container needs a unique identifier."; + throw GException::invalid_value(G_APPEND, msg); + } + + // Clone observation + GCTAOnOffObservation* ptr = obs.clone(); + + // Append to list + m_obs.push_back(ptr); + + // Return pointer + return ptr; +} + + +/***********************************************************************//** + * @brief Insert observation into container + * + * @param[in] index Observation index [0,...,size()-1]. + * @param[in] obs Observation. + * @return Pointer to deep copy of observation. + * + * @exception GException::out_of_range + * Observation index is out of range. + * @exception GException::invalid_value + * Observation with same instrument and identifier already + * exists in container. + * + * Inserts a deep copy of an observation into the container before the + * observation with the specified @p index. + ***************************************************************************/ +GCTAOnOffObservation* GCTAOnOffObservations::insert(const int& index, const GCTAOnOffObservation& obs) +{ + // Compile option: raise exception if index is out of range + #if defined(G_RANGE_CHECK) + if (isempty()) { + if (index > 0) { + throw GException::out_of_range(G_INSERT, index, 0, size()-1); + } + } + else { + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_INSERT, index, 0, size()-1); + } + } + #endif + + // Raise an exception if an observation with specified instrument and + // identifier already exists + int inx = get_index(obs.instrument(), obs.id()); + if (inx != -1 && inx != index) { + std::string msg = + "Attempt to insert \""+obs.instrument()+"\" observation with" + " identifier \""+obs.id()+"\" in observation container before" + " index "+gammalib::str(index)+", but an observation with the" + " same attributes exists already at index "+gammalib::str(inx)+ + " in the container.\n" + "Every observation for a given instrument in the observation" + " container needs a unique identifier."; + throw GException::invalid_value(G_INSERT, msg); + } + + // Clone observation + GCTAOnOffObservation* ptr = obs.clone(); + + // Clone observation and insert into list + m_obs.insert(m_obs.begin()+index, ptr); + + // Return pointer + return ptr; +} + + +/***********************************************************************//** + * @brief Remove observation from container + * + * @param[in] index Observation index [0,...,size()-1]. + * + * @exception GException::out_of_range + * Observation index is out of range. + * + * Removes observation of specified @p index from the container. + ***************************************************************************/ +void GCTAOnOffObservations::remove(const int& index) +{ + // Compile option: If index is outside boundary then raise exception + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_REMOVE, index, 0, size()-1); + } + #endif + + // Erase observation component from container + m_obs.erase(m_obs.begin() + index); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Append observations from observation container + * + * @param[in] obs Observations. + * + * @exception GException::invalid_value + * Observation with same instrument and identifier already + * exists in container. + * + * Appends deep copies of observations to the container. + ***************************************************************************/ +void GCTAOnOffObservations::extend(const GCTAOnOffObservations& obs) +{ + // Do nothing if observation container is empty + if (!obs.isempty()) { + + // Get size. Note that we extract the size first to avoid an + // endless loop that arises when a container is appended to + // itself. + int num = obs.size(); + + // Reserve enough space + reserve(size() + num); + + // Loop over all observations, clone them and append them to the + // list + for (int i = 0; i < num; ++i) { + + // Raise an exception if an observation with specified + // instrument and identifier already exists + int inx = get_index(obs[i]->instrument(), obs[i]->id()); + if (inx != -1) { + std::string msg = + "Attempt to append \""+obs[i]->instrument()+"\"" + " observation with identifier \""+obs[i]->id()+"\"" + " to observation container, but an observation with the" + " same attributes exists already at index "+ + gammalib::str(inx)+" in the container.\n" + "Every observation for a given instrument in the" + " observation container needs a unique identifier."; + throw GException::invalid_value(G_EXTEND, msg); + } + + // Append observation to container + m_obs.push_back(obs[i]->clone()); + } + + } // endif: observation container was not empty + + // Return + return; +} + + +/***********************************************************************//** + * @brief Signals if observation exists + * + * @param[in] instrument Instrument. + * @param[in] id Observation identifier. + * @return True if observation with specified @p instrument and identifier + * @p id exists. + * + * Searches all observations for a match with the specified @p instrument + * and @p identifier. If the specified attributes have been found, true is + * returned. + ***************************************************************************/ +bool GCTAOnOffObservations::contains(const std::string& instrument, + const std::string& id) const +{ + // Get observation index + int index = get_index(instrument, id); + + // Return + return (index != -1); +} + + +/***********************************************************************//** + * @brief Load observations from XML file + * + * @param[in] filename XML filename. + * + * Loads observation from a XML file into the container. Please refer to the + * read() method for more information about the structure of the XML file. + ***************************************************************************/ +void GCTAOnOffObservations::load(const std::string& filename) +{ + // Clear any existing observations + clear(); + + // Load XML document + GXml xml(filename); + + // Read observations from XML document + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save observations into XML file + * + * @param[in] filename XML filename. + * + * Saves observations into a XML file. Please refer to the read() method for + * more information about the structure of the XML file. + ***************************************************************************/ +void GCTAOnOffObservations::save(const std::string& filename) const +{ + // Declare empty XML document + GXml xml; + + // Write observations into XML file + write(xml); + + // Save XML document + xml.save(filename); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Read observations from XML document + * + * @param[in] xml XML document. + * + * @exception GException::invalid_instrument + * Invalid instrument encountered in XML file. + * + * Reads observations from the first observation list that is found in the + * XML document. The decoding of the instrument specific observation + * definition is done within the observation's GCTAOnOffObservation::read() method. + * The following file structure is expected: + * + * + * + * ... + * + * + * ... + * + * ... + * + * + * The @p name and @p id attributes allow for a unique identification of an + * observation within the observation container. The @p instrument + * attributes specifies the instrument to which the observation applies. + * At the moment there is only CTA but more instruments may be implemented + * later and require the need for a registry (see GObservations.cpp). + * + * The structure within the @p observation tag is defined by the instrument + * specific GCTAOnOffObservation class. + * + ***************************************************************************/ +void GCTAOnOffObservations::read(const GXml& xml) +{ + // Get pointer on observation library + const GXmlElement* lib = xml.element("observation_list", 0); + + // Loop over all observations + int n = lib->elements("observation"); + for (int i = 0; i < n; ++i) { + + // Get pointer on observation + const GXmlElement* obs = lib->element("observation", i); + + // Get attributes + std::string name = obs->attribute("name"); + std::string id = obs->attribute("id"); + std::string instrument = obs->attribute("instrument"); + + // Allocate observation + // (only CTA at the moment, implement registry if more) + GCTAOnOffObservation* ptr = new GCTAOnOffObservation; + + // If observation is valid then read its definition from XML file + if (ptr != NULL) { + + // Read definition + ptr->read(*obs); + + // Set attributes + ptr->name(name); + ptr->id(id); + + } // endif: observation was valid + + // ... otherwise throw an exception + else { + std::string msg = + "Problem allocating GCTAOnOffObservation object."; + throw GException::invalid_instrument(G_READ, msg); + } + + // Append observation to container + append(*ptr); + + // Free observation (the append method made a deep copy) + delete ptr; + + } // endfor: looped over all observations + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write observations into XML document + * + * @param[in] xml XML document. + * + * Write observations into the first observation library that is found in the + * XML document. In case that no observation library exists, one is added to + * the document. Please refer to the read() method for more information + * about the structure of the XML document. + ***************************************************************************/ +void GCTAOnOffObservations::write(GXml& xml) const +{ + // If there is no observation library then append one + if (xml.elements("observation_list") == 0) { + xml.append(GXmlElement("observation_list title=\"observation list\"")); + } + + // Get pointer on observation library + GXmlElement* lib = xml.element("observation_list", 0); + + // Write all observations into library + for (int i = 0; i < size(); ++i) { + + // Initialise pointer on observation + GXmlElement* obs = NULL; + + // Search corresponding observation + int n = xml.elements("observation"); + for (int k = 0; k < n; ++k) { + GXmlElement* element = xml.element("observation", k); + if (element->attribute("name") == m_obs[i]->name() && + element->attribute("id") == m_obs[i]->id() && + element->attribute("instrument") == m_obs[i]->instrument()) { + obs = element; + break; + } + } + + // If no observation with corresponding name, ID and instrument was + // found then append one now + if (obs == NULL) { + obs = lib->append("observation"); + obs->attribute("name", m_obs[i]->name()); + obs->attribute("id", m_obs[i]->id()); + obs->attribute("instrument", m_obs[i]->instrument()); + } + + // Write now observation + m_obs[i]->write(*obs); + + } // endfor: looped over all observaitons + + // Return + return; +} + + +/***********************************************************************//** + * @brief Load models from XML file + * + * @param[in] filename XML filename. + * + * Loads the models from a XML file. Please refer to the GModels::read() + * method for more information about the expected structure of the XML + * file. + ***************************************************************************/ +void GCTAOnOffObservations::models(const std::string& filename) +{ + // Load models + m_models.load(filename); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print observation list information + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing observation list information + ***************************************************************************/ +std::string GCTAOnOffObservations::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GCTAOnOffObservations ==="); + + // Append information + result.append("\n"+gammalib::parformat("Number of observations")); + result.append(gammalib::str(size())); + + // NORMAL: Append observations + if (chatter >= NORMAL) { + for (int i = 0; i < size(); ++i) { + result.append("\n"); + result.append((*this)[i]->print()); + } + } // endif: chatter was normal + + // EXPLICIT: Append models + if (chatter >= EXPLICIT) { + result.append("\n"+m_models.print()); + } // endif: chatter was explicit + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GCTAOnOffObservations::init_members(void) +{ + // Initialise members + m_obs.clear(); + m_models.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] obs Observation container. + * + * Copies all class members. Deep copies are created for the observations, + * which is what is needed to have a real physical copy of the members. + ***************************************************************************/ +void GCTAOnOffObservations::copy_members(const GCTAOnOffObservations& obs) +{ + // Copy attributes + m_models = obs.m_models; + + // Copy observations + m_obs.clear(); + for (int i = 0; i < obs.m_obs.size(); ++i) { + m_obs.push_back((obs.m_obs[i]->clone())); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + * + * Deallocates all observations. Since container classes that hold pointers + * need to handle the proper deallocation of memory, we loop here over all + * pointers and make sure that we deallocate all observations in the + * container. + ***************************************************************************/ +void GCTAOnOffObservations::free_members(void) +{ + // Free observations + for (int i = 0; i < m_obs.size(); ++i) { + delete m_obs[i]; + m_obs[i] = NULL; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Return observation index by instrument and identifier + * + * @param[in] instrument Instrument. + * @param[in] id Observation identifier. + * @return Observation index (-1 of instrument and observation identifier + * has not been found) + * + * Returns observation index based on the specified @p instrument and + * observation identifier @p id. If no observation with the specified + * attributes has been found, the method returns -1. + ***************************************************************************/ +int GCTAOnOffObservations::get_index(const std::string& instrument, + const std::string& id) const +{ + // Initialise index + int index = -1; + + // Search observation with specified instrument and id + for (int i = 0; i < size(); ++i) { + if ((m_obs[i]->instrument() == instrument) && + (m_obs[i]->id() == id)) { + index = i; + break; + } + } + + // Return index + return index; +} diff --git a/pyext/GArf.i b/pyext/GArf.i new file mode 100644 index 0000000..d3ce275 --- /dev/null +++ b/pyext/GArf.i @@ -0,0 +1,87 @@ +/*************************************************************************** + * GArf.i - XSPEC Auxiliary Response File class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GArf.i + * @brief XSPEC Auxiliary Response File class definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GArf.hpp" +#include "GTools.hpp" +%} + + +/***********************************************************************//** + * @class GArf + * + * @brief Auxiliary Response File class + ***************************************************************************/ +class GArf : public GBase { +public: + // Constructors and destructors + GArf(void); + explicit GArf(const std::string& filename); + explicit GArf(const GEbounds& ebds); + GArf(const GArf& arf); + virtual ~GArf(void); + + // Methods + void clear(void); + GArf* clone(void) const; + int size(void) const; + 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 write(GFits& fits) const; +}; + + +/***********************************************************************//** + * @brief GArf class extension + ***************************************************************************/ +%extend GArf { + double __getitem__(const int& index) { + if (index >= 0 && index < self->size()) { + return (*self)[index]; + } + else { + throw GException::out_of_range("__getitem__(int)", index, self->size()); + } + } + void __setitem__(const int& index, const double& value) { + if (index>=0 && index < self->size()) { + (*self)[index] = value; + return; + } + else { + throw GException::out_of_range("__setitem__(int)", index, self->size()); + } + } + int __len__() { + return (self->size()); + } + GArf copy() { + return (*self); + } +}; diff --git a/pyext/GPha.i b/pyext/GPha.i new file mode 100644 index 0000000..5d41719 --- /dev/null +++ b/pyext/GPha.i @@ -0,0 +1,92 @@ +/*************************************************************************** + * GPha.i - XSPEC Pulse Height Analyzer class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GPha.i + * @brief XSPEC Pulse Height Analyzer class + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GPha.hpp" +#include "GTools.hpp" +%} + + +/***********************************************************************//** + * @class GPha + * + * @brief Pulse Height Analyzer class + ***************************************************************************/ +class GPha : public GBase { +public: + // Constructors and destructors + GPha(void); + explicit GPha(const std::string& filename); + explicit GPha(const GEbounds& ebds); + GPha(const GPha& pha); + virtual ~GPha(void); + + // Methods + void clear(void); + GPha* clone(void) const; + int size(void) const; + double& at(const int& index); + const GEbounds& ebounds(void) const; + double counts(void) const; + const double& underflow(void) const; + const double& overflow(void) const; + 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 write(GFits& fits) const; +}; + + +/***********************************************************************//** + * @brief GPha class extension + ***************************************************************************/ +%extend GPha { + double __getitem__(const int& index) { + if (index >= 0 && index < self->size()) { + return (*self)[index]; + } + else { + throw GException::out_of_range("__getitem__(int)", index, self->size()); + } + } + void __setitem__(const int& index, const double& value) { + if (index>=0 && index < self->size()) { + (*self)[index] = value; + return; + } + else { + throw GException::out_of_range("__setitem__(int)", index, self->size()); + } + } + int __len__() { + return (self->size()); + } + GPha copy() { + return (*self); + } +}; diff --git a/pyext/GRmf.i b/pyext/GRmf.i new file mode 100644 index 0000000..9524e88 --- /dev/null +++ b/pyext/GRmf.i @@ -0,0 +1,82 @@ +/*************************************************************************** + * GRmf.i - XSPEC Redistribution Matrix File class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GRmf.i + * @brief XSPEC Redistribution Matrix File class + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GRmf.hpp" +#include "GTools.hpp" +%} + +/* __ Includes ___________________________________________________________ */ +%include "GTypemaps.i" + + +/***********************************************************************//** + * @class GRmf + * + * @brief Redistribution Matrix File class + ***************************************************************************/ +class GRmf : public GBase { +public: + // Constructors and destructors + GRmf(void); + explicit GRmf(const std::string& filename); + explicit GRmf(const GEbounds& etrue, const GEbounds& emeasured); + GRmf(const GRmf& rmf); + virtual ~GRmf(void); + + // Methods + void clear(void); + GRmf* clone(void) const; + int size(void) const; + int ntrue(void) const; + int nmeasured(void) const; + double& at(const int& itrue, const int& imeasured); + 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 write(GFits& fits) const; +}; + + +/***********************************************************************//** + * @brief GRmf class extension + ***************************************************************************/ +%extend GRmf { + double __getitem__(int GTuple[2]) { + return self->at(GTuple[0], GTuple[1]); + } + void __setitem__(int GTuple[2], const double& value) { + self->at(GTuple[0], GTuple[1]) = value; + } + int __len__() { + return (self->size()); + } + GRmf copy() { + return (*self); + } +}; diff --git a/pyext/GSkyRegion.i b/pyext/GSkyRegion.i new file mode 100644 index 0000000..81678fa --- /dev/null +++ b/pyext/GSkyRegion.i @@ -0,0 +1,63 @@ +/*************************************************************************** + * GSkyRegion.i - Sky region class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GSkyRegion.i + * @brief Sky region class SWIG file. + * @author Michael Mayer + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GSkyRegion.hpp" +#include "GTools.hpp" +%} + + +/***********************************************************************//** + * @class GSkyRegion + * + * @brief GSkyRegion base abstract class + ***************************************************************************/ +class GSkyRegion : public GBase { +public: + // Constructors and destructors + GSkyRegion(void); + GSkyRegion(const GSkyRegion& region); + virtual ~GSkyRegion(void); + + // Pure virtual methods + virtual void clear(void) = 0; + virtual GSkyRegion* clone(void) const = 0; + virtual void read(const std::string regstring) = 0; + virtual std::string write(void) const = 0; + virtual bool contains(const GSkyDir& dir) const = 0; + virtual bool overlaps(const GSkyRegion& reg) const = 0; + virtual bool contains(const GSkyRegion& reg) const = 0; + + // Implemented methods + std::string type(void) const; + std::string name(void) const; + double solidangle(void) const; +}; + + +/***********************************************************************//** + * @brief GSkyRegion class extension + ***************************************************************************/ diff --git a/pyext/GSkyRegionCircle.i b/pyext/GSkyRegionCircle.i new file mode 100644 index 0000000..4b57312 --- /dev/null +++ b/pyext/GSkyRegionCircle.i @@ -0,0 +1,71 @@ +/*************************************************************************** + * GSkyRegionCircle.i - Sky region class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GSkyRegionCircle.i + * @brief Sky region class SWIG file. + * @author Michael Mayer + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GSkyRegionCircle.hpp" +#include "GTools.hpp" +%} + + +/***********************************************************************//** + * @class GSkyRegionCircle + * + * @brief GSkyRegionCircle dervied sky region class + ***************************************************************************/ +class GSkyRegionCircle : public GSkyRegion { +public: + // Constructors and destructors + GSkyRegionCircle(void); + GSkyRegionCircle(const GSkyRegionCircle& circle); + explicit GSkyRegionCircle(const std::string &line); + GSkyRegionCircle(GSkyDir ¢re, const double &radius); + // GSkyRegionCircle(const double &ra, const double &dec, const double &radius); + virtual ~GSkyRegionCircle(void); + + // Implemented methods + void clear(void); + GSkyRegionCircle* clone(void) const; + GSkyDir centre(void) const; + double radius(void) const; + void radius(const double& radius); + void centre(const GSkyDir& centre); + void centre(const double& ra,const double& dec); + void read(const std::string line); + std::string write(void) const; + bool contains(const GSkyDir& dir) const; + bool contains(const GSkyRegion& reg) const; + bool overlaps(const GSkyRegion& reg) const; +}; + + +/***********************************************************************//** + * @brief GSkyRegionCircle class extension + ***************************************************************************/ +%extend GSkyRegionCircle { + GSkyRegionCircle copy() { + return (*self); + } +}; diff --git a/pyext/GSkyRegions.i b/pyext/GSkyRegions.i new file mode 100644 index 0000000..56b72af --- /dev/null +++ b/pyext/GSkyRegions.i @@ -0,0 +1,95 @@ +/*************************************************************************** + * GSkyRegions.i - Sky region container class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GSkyRegions.i + * @brief Sky regions container class SWIG file + * @author Pierrick Martin + */ + +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GBase.hpp" +#include "GSkyRegion.hpp" +%} + +/***********************************************************************//** + * @class GSkyRegions + * + * @brief Sky region container class + ***************************************************************************/ +class GSkyRegions : public GBase { +public: + // Constructors and destructors + GSkyRegions(void); + GSkyRegions(const GSkyRegions& regions); + explicit GSkyRegions(const std::string& filename); + virtual ~GSkyRegions(void); + + // Methods + void clear(void); + GSkyRegions* clone(void) const; + GSkyRegion* at(const int& index); + const GSkyRegion* at(const int& index) const; + int size(void) const; + bool isempty(void) const; + GSkyRegion* set(const int& index, const GSkyRegion& region); + GSkyRegion* set(const std::string& name, const GSkyRegion& region); + GSkyRegion* append(const GSkyRegion& region); + GSkyRegion* insert(const int& index, const GSkyRegion& region); + GSkyRegion* insert(const std::string& name, const GSkyRegion& region); + void remove(const int& index); + void remove(const std::string& name); + void reserve(const int& num); + void extend(const GSkyRegions& regions); + bool contains(const std::string& name) const; + void load(const std::string& filename); + void save(const std::string& filename) const; + +}; + +/***********************************************************************//** + * @brief GSkyRegions class extension + ***************************************************************************/ +%extend GSkyRegions { + GSkyRegion* __getitem__(const int& index) { + if (index >= 0 && index < self->size()) { + return (*self)[index]; + } + else { + throw GException::out_of_range("__getitem__(int)", index, self->size()); + } + } + void __setitem__(const int& index, const GSkyRegion& val) { + if (index>=0 && index < self->size()) { + self->set(index, val); + return; + } + else { + throw GException::out_of_range("__setitem__(int)", index, self->size()); + } + } + int __len__() { + return (self->size()); + } + GSkyRegions copy() { + return (*self); + } +}; diff --git a/pyext/Makefile.am b/pyext/Makefile.am index 49247bf..fb087eb 100755 --- a/pyext/Makefile.am +++ b/pyext/Makefile.am @@ -62,6 +62,7 @@ BUILT_SOURCES = gammalib/app_wrap.cpp gammalib/app.py \ gammalib/sky_wrap.cpp gammalib/sky.py \ gammalib/support_wrap.cpp gammalib/support.py \ gammalib/xml_wrap.cpp gammalib/xml.py \ + gammalib/xspec_wrap.cpp gammalib/xspec.py \ gammalib/test_wrap.cpp gammalib/test.py \ $(INST_SWIG_TARGETS) diff --git a/pyext/gammalib/sky.i b/pyext/gammalib/sky.i index e2603ae..748f85d 100644 --- a/pyext/gammalib/sky.i +++ b/pyext/gammalib/sky.i @@ -50,6 +50,8 @@ %include "GSkyDir.i" %include "GSkyPixel.i" %include "GSkymap.i" +%include "GSkyRegion.i" +%include "GSkyRegionCircle.i" %include "GWcs.i" %include "GWcsRegistry.i" %include "GWcslib.i" diff --git a/pyext/gammalib/xspec.i b/pyext/gammalib/xspec.i new file mode 100644 index 0000000..ee05ceb --- /dev/null +++ b/pyext/gammalib/xspec.i @@ -0,0 +1,52 @@ +/*************************************************************************** + * xspec.i - XSPEC module * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + * ----------------------------------------------------------------------- * + * Usage: * + * swig -c++ -python -Wall xspec.i * + ***************************************************************************/ +/** + * @file xspec.i + * @brief XSPEC module + * @author Juergen Knoedlseder + */ +%module xspec +%feature("autodoc", "1"); + +/* __ Headers needed for compilation _____________________________________ */ +%{ +#include +#include "GException.hpp" +#include "GTools.hpp" +%} + +/* __ Include standard typemaps for vectors and strings __________________ */ +%include stl.i + +/* __ Include interface classes __________________________________________ */ +%import(module="gammalib.base") "GBase.i"; +//%import(module="gammalib.base") "GContainer.i"; + +/* __ Make sure that exceptions are catched ______________________________ */ +%import(module="gammalib.support") "GException.i"; + +/* __ Application module _________________________________________________ */ +%include "GArf.i" +%include "GPha.i" +%include "GRmf.i" diff --git a/pyext/setup.py.in b/pyext/setup.py.in index 2ff57e2..6819d3b 100644 --- a/pyext/setup.py.in +++ b/pyext/setup.py.in @@ -1,7 +1,7 @@ ############################################################################# # Disutils setup file for gammalib Python extension module # # ------------------------------------------------------------------------- # -# Copyright (C) 2011-2012 Juergen Knoedlseder # +# Copyright (C) 2011-2013 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -31,7 +31,7 @@ pkg_name = 'gammalib' # Define core modules core_modules = ['app', 'base', 'fits', 'linalg', 'model', 'numerics', \ - 'obs', 'opt', 'sky', 'support', 'test', 'xml'] + 'obs', 'opt', 'sky', 'support', 'test', 'xml', 'xspec'] # Detect execution mode is_build = True diff --git a/src/Makefile.am b/src/Makefile.am index 38e23d0..8dd9157 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -16,6 +16,7 @@ SUBDIRS = base \ fits \ xml \ vo \ + xspec \ sky \ opt \ obs \ @@ -59,6 +60,7 @@ libgamma_la_LIBADD = base/libbase.la \ fits/libfits.la \ xml/libxml.la \ vo/libvo.la \ + xspec/libxspec.la \ sky/libsky.la \ opt/libopt.la \ obs/libobs.la \ diff --git a/src/obs/GEbounds.cpp b/src/obs/GEbounds.cpp index e7b8235..ff3f32b 100644 --- a/src/obs/GEbounds.cpp +++ b/src/obs/GEbounds.cpp @@ -699,8 +699,14 @@ void GEbounds::write(GFits* file, const std::string& extname) const * @param[in] eng Energy. * @return Bin index. * - * Returns the energy boundary bin index for a given energy. If the energy - * falls outside all boundaries, -1 is returned. + * Returns the energy boundary bin index for a given energy. By convention, + * the limits for an energy bin are defined as + * + * min <= energy < max + * + * i.e. and energy equals to max falls above the largest energy. + * + * If the energy falls outside all boundaries, -1 is returned. ***************************************************************************/ int GEbounds::index(const GEnergy& eng) const { @@ -709,7 +715,7 @@ int GEbounds::index(const GEnergy& eng) const // Search all energy boundaries for containment for (int i = 0; i < m_num; ++i) { - if (eng >= m_min[i] && eng <= m_max[i]) { + if (eng >= m_min[i] && eng < m_max[i]) { index = i; break; } diff --git a/src/sky/GSkyRegion.cpp b/src/sky/GSkyRegion.cpp new file mode 100644 index 0000000..c51e664 --- /dev/null +++ b/src/sky/GSkyRegion.cpp @@ -0,0 +1,172 @@ +/*************************************************************************** + * GSkyRegion.cpp - Abstract virtual sky region base class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GSkyRegion.cpp + * @brief Abstract sky region base class interface definition + * @author Pierrick Martin + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include +#include "GSkyRegion.hpp" + +/* __ Method name definitions ____________________________________________ */ + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GSkyRegion::GSkyRegion(void) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] region sky region. + ***************************************************************************/ +GSkyRegion::GSkyRegion(const GSkyRegion& region) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(region); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GSkyRegion::~GSkyRegion(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] region sky region. + * @return sky region. + ***************************************************************************/ +GSkyRegion& GSkyRegion::operator=(const GSkyRegion& region) +{ + // Execute only if object is not identical + if (this != ®ion) { + + // Free members + free_members(); + + // Initialise private members for clean destruction + init_members(); + + // Copy members + copy_members(region); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GSkyRegion::init_members(void) +{ + // Initialise members + m_type.clear(); + m_name.clear(); + m_solid = 0.0; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] region sky region. + ***************************************************************************/ +void GSkyRegion::copy_members(const GSkyRegion& region) +{ + // Copy members + m_type = region.m_type; + m_name = region.m_name; + m_solid = region.m_solid; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GSkyRegion::free_members(void) +{ + // Return + return; +} diff --git a/src/sky/GSkyRegionCircle.cpp b/src/sky/GSkyRegionCircle.cpp new file mode 100644 index 0000000..f2d644a --- /dev/null +++ b/src/sky/GSkyRegionCircle.cpp @@ -0,0 +1,544 @@ +/*************************************************************************** + * GSkyRegionCircle.cpp - Circular sky region class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GSkyRegionCircle.hpp + * @brief Circular sky region implementation + * @author Michael Mayer + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GSkyRegionCircle.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_CONSTRUCTOR "GSkyRegionCircle::GSkyRegionCircle(GSkyDir&, double&)" +#define G_RADIUS "GSkyRegionCircle::radius(double&)" +#define G_READ "GSkyRegionCircle::read(std::string&)" +#define G_CONTAINS "GSkyRegionCircle::contains(GSkyRegion&)" +#define G_OVERLAPS "GSkyRegionCircle::overlaps(GSkyRegion&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + +/* __ Prototype __________________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GSkyRegionCircle::GSkyRegionCircle(void) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Direction constructor + * + * @param[in] centre GSkyDir of Circle center. + * @param[in] radius Radius of the region [deg]. + * + * @exception GException::invalid_argument + * Radius value is less than 0. + ***************************************************************************/ +GSkyRegionCircle::GSkyRegionCircle(GSkyDir& centre, const double& radius) +{ + // Initialise members + init_members(); + + // set members + m_centre = centre; + m_radius = radius; + + // Check if radius is valid + if (m_radius < 0.0) { + throw GException::invalid_argument(G_CONSTRUCTOR, + "Radius of a region can't be less than 0"); + } + + // compute solid angle + compute_solid(); + + // Return + return; +} + +/***********************************************************************//** + * @brief string constructor + * + * @param[in] line a string containing a line of a ds9 region file + ***************************************************************************/ +GSkyRegionCircle::GSkyRegionCircle(const std::string& line) +{ + // Initialise members + init_members(); + + // read from line + read(line); + + // Return + return; + +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] circle Circular sky region. + ***************************************************************************/ +GSkyRegionCircle::GSkyRegionCircle(const GSkyRegionCircle& circle) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(circle); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GSkyRegionCircle::~GSkyRegionCircle(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] circle Circular sky region. + * @return Circular sky region. + ***************************************************************************/ +GSkyRegionCircle& GSkyRegionCircle::operator= (const GSkyRegionCircle& circle) +{ + // Execute only if object is not identical + if (this != &circle) { + + // Free members + free_members(); + + // Initialise private members for clean destruction + init_members(); + + // Copy members + copy_members(circle); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear instance + ***************************************************************************/ +void GSkyRegionCircle::clear(void) +{ + // Free members + free_members(); + this->GSkyRegion::free_members(); + + // Initialise private members + this->GSkyRegion::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone object + ***************************************************************************/ +GSkyRegionCircle* GSkyRegionCircle::clone(void) const +{ + // Clone this Region + return new GSkyRegionCircle(*this); +} + + +/***********************************************************************//** + * @brief Set radius value of region + * + * @param[in] radius Radius value [deg]. + * + * @exception GException::invalid_argument + * Radius value is less than 0. + ***************************************************************************/ +void GSkyRegionCircle::radius(const double& radius) +{ + // Set radius value + m_radius = radius; + + if (m_radius < 0.0) { + throw GException::invalid_argument(G_RADIUS, + "Radius of a region can't be less than 0"); + } + + // Recompute the solid angle + compute_solid(); +} + + +/***********************************************************************//** + * @brief Read region from string + * + * @param[in] line String in ds9 format. + * + * @exception GException::invalid_value + * Invalid value found in ds9 format string. + ***************************************************************************/ +void GSkyRegionCircle::read(const std::string& line) +{ + // Clear the current instance + clear(); + + // Split the string into 2 parts seperated by # + std::vector substrings = gammalib::split(line,"#"); + std::string region_def = (substrings.size() > 0) ? substrings[0] : ""; + std::string comment = (substrings.size() > 1) ? substrings[1] : ""; + + // Finding the circle + if (region_def.find("circle") == std::string::npos) { + throw GException::invalid_value(G_READ, + "Cannnot find the key word \"circle\" in provided string"); + } + + // Get the the coordinate system of the values + std::string system = gammalib::split(region_def,";")[0]; + + // Get the substring of the important values + unsigned pos = region_def.find("circle("); + unsigned end = region_def.find(")"); + std::string circlestring = region_def.substr(pos+7,end); + circlestring.erase(circlestring.find(")"),1); + + // Get the values of the region x,y,and radius + std::vector values = gammalib::split(circlestring,","); + if (values.size()!=3) { + throw GException::invalid_value(G_READ, + "Invalid number of arguments, circle takes exactly 3 arguments"); + } + double x = gammalib::todouble(values[0]); + double y = gammalib::todouble(values[1]); + double radius = gammalib::todouble(values[2]); + + // Initialise sky dir + GSkyDir dir = GSkyDir(); + + // Set the values in the correct system + if (system == "fk5") { + dir.radec_deg(x,y); + } + else if (system == "galactic") { + dir.lb_deg(x,y); + } + else { + std::stringstream s; + s << "Provided coordinate system \"" << system << "\" is not supported."; + throw GException::invalid_value(G_READ, s.str()); + } + + // Setting par values + m_centre = dir; + m_radius = radius; + + // Throwing exception if radius is less than 0 + if (m_radius < 0.0) { + throw GException::invalid_value(G_READ, + "Radius of a region can't be less than 0"); + } + + // Compute solid angle + compute_solid(); + + // Check if there is a given name for the region and set it + std::vectorcomments = gammalib::split(comment, " "); + for (int i = 0; i < comments.size(); i++) { + if (gammalib::contains(comments[i], "text")) { + std::vector attributes = gammalib::split(comments[i], "="); + if (attributes.size() < 2) { + throw GException::invalid_value(G_READ, + "Invalid number of arguments, type of attribute must" + " be key=value"); + } + m_name = attributes[1]; + } + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief write region into a string + * + * @return String to be written in a ds9 region file + ***************************************************************************/ +std::string GSkyRegionCircle::write(void) const +{ + std::string result; + result.append("fk5;circle("); + result.append(gammalib::str(m_centre.ra_deg())); + result.append(","); + result.append(gammalib::str(m_centre.dec_deg())); + result.append(","); + result.append(gammalib::str(m_radius)); + result.append(") # text="); + result.append(m_name); + result.append("\n"); + return result; +} + + +/***********************************************************************//** + * @brief Print circular region + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing region information. + ***************************************************************************/ +std::string GSkyRegionCircle::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result;; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append string + result.append("=== GSkyRegionCircle ==="); + result.append("\n("); + result.append(gammalib::str(m_centre.ra_deg())); + result.append(","); + result.append(gammalib::str(m_centre.dec_deg())); + result.append(","); + result.append(gammalib::str(m_radius)); + result.append(")"); + + } // endif: chatter was not silent + + // Return result + return result; +} + +/***********************************************************************//** + * @brief checks if coordinate is within region + * + * @param[in] dir Sky direction. + ***************************************************************************/ +bool GSkyRegionCircle::contains(const GSkyDir& dir) const +{ + // Initialise return value + bool dir_is_in = false; + + // calculate distance from dir to centre + double distance = dir.dist_deg(m_centre); + + // Check if distance < radius + if (distance <= radius()) { + dir_is_in = true; + } + + // Return bool + return dir_is_in; +} + +/***********************************************************************//** + * @brief checks if region is fully contained within this region + * + * @param[in] reg Sky region. + * + * @exception GException::feature_not_implemented + * Regions differ in type. + ***************************************************************************/ +bool GSkyRegionCircle::contains(const GSkyRegion& reg) const +{ + // Initialise return value + bool fully_inside = false; + + // If other region is circle use a simple way to calculate + if (reg.type() == "Circle") { + + // create circular region from reg + GSkyRegionCircle* regcirc = (GSkyRegionCircle*)reg.clone(); + + // calculate angular distance between the centers + double ang_dist = m_centre.dist_deg(regcirc->centre()); + + // Check if the region is contained in this + if ((ang_dist + regcirc->radius()) <= m_radius) { + fully_inside = true; + } + + } + + // ... otherwise throw an exception + else { + throw GException::feature_not_implemented(G_CONTAINS, + "Cannot compare two different region types yet"); + } + + // Return value + return fully_inside; +} + +/***********************************************************************//** + * @brief Checks if region is overlapping with this region + * + * @param[in] reg Sky region. + * + * @exception GException::feature_not_implemented + * Regions differ in type. + ***************************************************************************/ +bool GSkyRegionCircle::overlaps(const GSkyRegion& reg) const +{ + // Initialise return value + bool overlap = false; + + // If other region is circle use a simple way to calculate + if (reg.type() == "Circle") { + + // Create circular region from reg + GSkyRegionCircle* regcirc = (GSkyRegionCircle*)reg.clone(); + + // Calculate angular distance between the centers + double ang_dist = m_centre.dist_deg(regcirc->centre()); + + // Check if the distance is smaller than the sum of both radii + if (ang_dist <= (m_radius + regcirc->radius())) { + overlap = true; + } + + } + + // ... otherwise throw an exception + else { + throw GException::feature_not_implemented(G_OVERLAPS, + "Cannot compare two different region types yet"); + } + + // Return value + return overlap; +} + + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GSkyRegionCircle::init_members(void) +{ + // Initialise members + m_centre = GSkyDir(); + m_radius = 0.0; + m_type = "Circle"; + + //Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] region Circular sky region. + ***************************************************************************/ +void GSkyRegionCircle::copy_members(const GSkyRegionCircle& region) +{ + // Copy attributes + m_centre = region.m_centre; + m_radius = region.m_radius; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GSkyRegionCircle::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief compute solid angle + ***************************************************************************/ +void GSkyRegionCircle::compute_solid(void) +{ + // convert radius to radians + double radius_rad = m_radius * gammalib::deg2rad; + + // compute solid angle + m_solid = gammalib::twopi * (1.0 - std::cos(radius_rad)); + + // Return + return; +} diff --git a/src/sky/GSkyRegions.cpp b/src/sky/GSkyRegions.cpp new file mode 100644 index 0000000..23eb4b4 --- /dev/null +++ b/src/sky/GSkyRegions.cpp @@ -0,0 +1,950 @@ +/*************************************************************************** + * GSkyRegions.cpp - Sky region container class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Pierrick Martin * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** +* @file GSkyRegions.cpp +* @brief Sky regions container class definition +* @author Pierrick Martin +*/ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#include +#endif +#include "GBase.hpp" +#include "GSkyRegion.hpp" +#include "GSkyRegionCircle.hpp" +#include "GSkyRegions.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_ACCESS "GSkyRegions::operator[](std::string&)" +#define G_AT "GSkyRegions::at(int&)" +#define G_SET1 "GSkyRegions::set(int&, GSkyRegion&)" +#define G_SET2 "GSkyRegions::set(std::string&, GSkyRegion&)" +#define G_APPEND "GSkyRegions::append(GSkyRegion&)" +#define G_INSERT1 "GSkyRegions::insert(int&, GSkyRegion&)" +#define G_INSERT2 "GSkyRegions::insert(std::string&, GSkyRegion&)" +#define G_REMOVE1 "GSkyRegions::remove(int&)" +#define G_REMOVE2 "GSkyRegions::remove(std::string&)" +#define G_EXTEND "GSkyRegions::extend(GSkyRegions&)" +#define G_LOAD "GSkyRegions::load(std::string&)" +#define G_SAVE "GSkyRegions::save(std::string&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GSkyRegions::GSkyRegions(void) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] regions region container. + ***************************************************************************/ +GSkyRegions::GSkyRegions(const GSkyRegions& regions) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(regions); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Load constructor + * + * @param[in] filename DS9 region file. + * + * Constructs region container from a DS9 region file. + ***************************************************************************/ +GSkyRegions::GSkyRegions(const std::string& filename) +{ + // Initialise members + init_members(); + + // Load XML file + load(filename); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GSkyRegions::~GSkyRegions(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] regions region container. + * @return region container. + ***************************************************************************/ +GSkyRegions& GSkyRegions::operator=(const GSkyRegions& regions) +{ + // Execute only if object is not identical + if (this != ®ions) { + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(regions); + + } // endif: object was not identical + + // Return + return *this; +} + + +/***********************************************************************//** + * @brief Return pointer to region + * + * @param[in] name region name. + * + * @exception GException::invalid_argument + * region with specified name not found in container. + * + * Returns a pointer to the region with the specified @p name. + ***************************************************************************/ +GSkyRegion* GSkyRegions::operator[](const std::string& name) +{ + // Get region index + int index = get_index(name); + + // Throw exception if region name was not found + if (index == -1) { + std::string msg ="Region with name \""+name+"\" not found in container."; + throw GException::invalid_argument(G_ACCESS, msg); + } + + // Return pointer + return m_regions[index]; +} + + +/***********************************************************************//** + * @brief Return pointer to region (const version) + * + * @param[in] name region name. + * + * @exception GException::invalid_argument + * region with specified name not found in container. + * + * Returns a const pointer to the region with the specified @p name. + ***************************************************************************/ +const GSkyRegion* GSkyRegions::operator[](const std::string& name) const +{ + // Get region index + int index = get_index(name); + + // Throw exception if region name was not found + if (index == -1) { + std::string msg ="Region with name \""+name+"\" not found in container."; + throw GException::invalid_argument(G_ACCESS, msg); + } + + // Return pointer + return m_regions[index]; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear object + * + * Removes all regions from the container. + ***************************************************************************/ +void GSkyRegions::clear(void) +{ + // Free class members + free_members(); + + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone instance + * + * @return Pointer to deep copy of region container + * + * Makes a deep copy of the region container instance. + ***************************************************************************/ +GSkyRegions* GSkyRegions::clone(void) const +{ + return new GSkyRegions(*this); +} + + +/***********************************************************************//** + * @brief Return pointer to region + * + * @param[in] index region index [0,...,size()-1]. + * + * @exception GException::out_of_range + * region index is out of range. + * + * Returns a pointer to the region with the specified @p index. + ***************************************************************************/ +GSkyRegion* GSkyRegions::at(const int& index) +{ + // Raise exception if index is out of range + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_AT, index, 0, size()-1); + } + + // Return pointer + return m_regions[index]; +} + + +/***********************************************************************//** + * @brief Return pointer to region (const version) + * + * @param[in] index region index [0,...,size()-1]. + * + * @exception GException::out_of_range + * region index is out of range. + * + * Returns a const pointer to the region with the specified @p index. + ***************************************************************************/ +const GSkyRegion* GSkyRegions::at(const int& index) const +{ + // Raise exception if index is out of range + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_AT, index, 0, size()-1); + } + + // Return pointer + return m_regions[index]; +} + + +/***********************************************************************//** + * @brief Set region in container + * + * @param[in] index sky region index [0,...,size()-1]. + * @param[in] region sky region. + * @return Pointer to deep copy of sky region. + * + * @exception GException::out_of_range + * region index is out of range. + * @exception GException::invalid_value + * Name of region exists already in container. + * + * Set sky region in the container. A deep copy of the region will be made. + ***************************************************************************/ +GSkyRegion* GSkyRegions::set(const int& index, const GSkyRegion& region) +{ + // Compile option: raise exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_SET1, index, 0, size()-1); + } + #endif + + // Check if a region with specified name does not yet exist + int inx = get_index(region.name()); + if (inx != -1 && inx != index) { + std::string msg = + "Attempt to set region with name \""+region.name()+"\" in region" + " container at index "+gammalib::str(index)+", but a region with" + " the same name exists already at index "+gammalib::str(inx)+ + " in the container.\n" + "Every region in the region container needs a unique name."; + throw GException::invalid_value(G_SET1, msg); + } + + // Delete any existing region + if (m_regions[index] != NULL) delete m_regions[index]; + + // Assign new region by cloning + m_regions[index] = region.clone(); + + // Return pointer to region + return m_regions[index]; +} + + +/***********************************************************************//** + * @brief Set region in container + * + * @param[in] name region name. + * @param[in] region region pointer. + * @return Pointer to deep copy of region. + * + * @exception GException::invalid_argument + * region with specified name not found in container. + * @exception GException::invalid_value + * Name of region exists already in container. + * + * Set region in the container. A deep copy of the region will be made. + ***************************************************************************/ +GSkyRegion* GSkyRegions::set(const std::string& name, const GSkyRegion& region) +{ + // Get parameter index + int index = get_index(name); + + // Throw exception if parameter name was not found + if (index == -1) { + throw GException::invalid_argument(G_SET2, name); + } + + // Check if a region with specified name does not yet exist + int inx = get_index(region.name()); + if (inx != -1 && inx != index) { + std::string msg = + "Attempt to set region with name \""+region.name()+"\" in region" + " container at index "+gammalib::str(index)+", but a region with" + " the same name exists already at index "+gammalib::str(inx)+ + " in the container.\n" + "Every region in the region container needs a unique name."; + throw GException::invalid_value(G_SET2, msg); + } + + // Delete any existing region + if (m_regions[index] != NULL) delete m_regions[index]; + + // Assign new region by cloning + m_regions[index] = region.clone(); + + // Return pointer to region + return m_regions[index]; +} + + +/***********************************************************************//** + * @brief Append region to container + * + * @param[in] region region. + * @return Pointer to deep copy of region. + * + * @exception GException::invalid_value + * Name of region exists already in container. + * + * Appends region to the container by making a deep copy of the region and + * storing its pointer. + ***************************************************************************/ +GSkyRegion* GSkyRegions::append(const GSkyRegion& region) +{ + // Check if a region with specified name does not yet exist + int inx = get_index(region.name()); + if (inx != -1) { + std::string msg = + "Attempt to append region with name \""+region.name()+"\" to region" + " container, but a region with the same name exists already at" + " index "+gammalib::str(inx)+" in the container.\n" + "Every region in the region container needs a unique name."; + throw GException::invalid_value(G_APPEND, msg); + } + + // Create deep copy of region + GSkyRegion* ptr = region.clone(); + + // Append deep copy of region + m_regions.push_back(ptr); + + // Return pointer to region + return ptr; +} + + +/***********************************************************************//** + * @brief Insert region into container + * + * @param[in] index region index [0,...,size()-1]. + * @param[in] region region. + * @return Pointer to deep copy of region. + * + * @exception GException::out_of_range + * region index is out of range. + * @exception GException::invalid_value + * Name of region exists already in container. + * + * Inserts a @p region into the container before the region with the specified + * @p index. + ***************************************************************************/ +GSkyRegion* GSkyRegions::insert(const int& index, const GSkyRegion& region) +{ + // Compile option: raise exception if index is out of range + #if defined(G_RANGE_CHECK) + if (isempty()) { + if (index > 0) { + throw GException::out_of_range(G_INSERT1, index, 0, size()-1); + } + } + else { + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_INSERT1, index, 0, size()-1); + } + } + #endif + + // Check if a region with specified name does not yet exist + int inx = get_index(region.name()); + if (inx != -1) { + std::string msg = + "Attempt to insert region with name \""+region.name()+"\" in region" + " container before index "+gammalib::str(index)+", but a region" + " with the same name exists already at index "+gammalib::str(inx)+ + " in the container.\n" + "Every region in the region container needs a unique name."; + throw GException::invalid_value(G_INSERT1, msg); + } + + // Create deep copy of region + GSkyRegion* ptr = region.clone(); + + // Inserts deep copy of region + m_regions.insert(m_regions.begin()+index, ptr); + + // Return pointer to region + return ptr; +} + + +/***********************************************************************//** + * @brief Insert region into container + * + * @param[in] name region name. + * @param[in] region region. + * @return Pointer to deep copy of region. + * + * @exception GException::invalid_argument + * region with specified name not found in container. + * @exception GException::invalid_value + * Name of region exists already in container. + * + * Inserts a @p region into the container before the region with the specified + * @p name. + ***************************************************************************/ +GSkyRegion* GSkyRegions::insert(const std::string& name, const GSkyRegion& region) +{ + // Get parameter index + int index = get_index(name); + + // Throw exception if parameter name was not found + if (index == -1) { + throw GException::invalid_argument(G_INSERT2, name); + } + + // Check if a region with specified name does not yet exist + int inx = get_index(region.name()); + if (inx != -1) { + std::string msg = + "Attempt to insert region with name \""+region.name()+"\" in region" + " container before index "+gammalib::str(index)+", but a region" + " with the same name exists already at index "+gammalib::str(inx)+ + " in the container.\n" + "Every region in the region container needs a unique name."; + throw GException::invalid_value(G_INSERT2, msg); + } + + // Create deep copy of region + GSkyRegion* ptr = region.clone(); + + // Inserts deep copy of region + m_regions.insert(m_regions.begin()+index, ptr); + + // Return pointer to region + return ptr; +} + + +/***********************************************************************//** + * @brief Remove region from container + * + * @param[in] index region index [0,...,size()-1]. + * + * @exception GException::out_of_range + * region index is out of range. + * + * Remove region of specified @p index from container. + ***************************************************************************/ +void GSkyRegions::remove(const int& index) +{ + // Compile option: raise exception if index is out of range + #if defined(G_RANGE_CHECK) + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_REMOVE1, index, 0, size()-1); + } + #endif + + // Erase region component from container + m_regions.erase(m_regions.begin() + index); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Remove region from container + * + * @param[in] name region name. + * + * @exception GException::invalid_argument + * region with specified name not found in container. + * + * Remove region of specified @p name from container. + ***************************************************************************/ +void GSkyRegions::remove(const std::string& name) +{ + // Get parameter index + int index = get_index(name); + + // Throw exception if region name was not found + if (index == -1) { + throw GException::invalid_argument(G_REMOVE2, name); + } + + // Erase region component from container + m_regions.erase(m_regions.begin() + index); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Append region container + * + * @param[in] regions region container. + * + * Append region container to the container. + ***************************************************************************/ +void GSkyRegions::extend(const GSkyRegions& regions) +{ + // Do nothing if region container is empty + if (!regions.isempty()) { + + // Get size. Note that we extract the size first to avoid an + // endless loop that arises when a container is appended to + // itself. + int num = regions.size(); + + // Reserve enough space + reserve(size() + num); + + // Loop over all region components and append pointers to deep copies + for (int i = 0; i < num; ++i) { + + // Check if region name does not yet exist + int inx = get_index(regions[i]->name()); + if (inx != -1) { + std::string msg = + "Attempt to append region with name \""+regions[i]->name()+ + "\" to region container, but a region with the same name" + " exists already at index "+gammalib::str(inx)+" in the" + " container.\n" + "Every region in the region container needs a unique name."; + throw GException::invalid_value(G_EXTEND, msg); + } + + // Append region to container + m_regions.push_back(regions[i]->clone()); + + } // endfor: looped over all regions + + } // endif: region container was not empty + + // Return + return; +} + + +/***********************************************************************//** + * @brief Signals if region name exists + * + * @param[in] name region name. + * @return True if region with specified @p name exists. + * + * Searches all region names for a match with the specified @p name. If the + * specified name has been found, true is returned. + ***************************************************************************/ +bool GSkyRegions::contains(const std::string& name) const +{ + // Get region index + int index = get_index(name); + + // Return + return (index != -1); +} + + +/***********************************************************************//** + * @brief Load regions from DS9 region file + * + * @param[in] filename DS9 region filename. + * + * @exception GException::file_open_error + * File could not be opened. + * + * Loads all regions from a DS9 region file. See the read() method for more + * information about the expected structure of the XML file. + ***************************************************************************/ +void GSkyRegions::load(const std::string& filename) +{ + // Clear any existing regions + clear(); + + // Open file. Throw an exception if opening failed. + std::ifstream ds9file; + ds9file.open(filename.c_str()); + if (ds9file.is_open()) { + + // Loop over file lines + std::string fileline=""; + std::string coordsys="galactic"; + while ( ds9file.good() ) { + + // Read one line + getline (ds9file,fileline); + + // If line is a comment then continue + if (fileline[0] == '#') { + continue; + } + + // Check for global definition of coordinate system + if (std::string::npos != fileline.find("fk5")) { + coordsys="fk5"; + } + + // If region is a circle + if (std::string::npos != fileline.find("circle")) { + + // Create instance of GSkyRegion object + GSkyRegionCircle region; + + // If coordinate system and region defined on the same line + if ((std::string::npos != fileline.find("fk5")) || + (std::string::npos != fileline.find("galactic"))) { + region.read(fileline); + append(region); + } + + // else, prepend the coordinate system + else { + std::string newfileline=coordsys; + newfileline.append("; "); + newfileline.append(fileline); + region.read(newfileline); + append(region); + } + } + } + + // Close file + ds9file.close(); + + // Store filename + m_filename = filename; + + } + + // File could not be opened + else { + throw GException::file_open_error(G_LOAD, filename); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save regions into DS9 region file + * + * @param[in] filename DS9 region filename. + * + * @exception GException::file_open_error + * File could not be opened. + * + * Saves all regions in the container into a DS9 region file. + ***************************************************************************/ +void GSkyRegions::save(const std::string& filename) const +{ + // Open file + std::ofstream ds9file; + ds9file.open(filename.c_str()); + + // If file opened correctly, then save regions + if (ds9file.is_open()) { + + // Write global definition + std::string fileline; + fileline.append("# Region file format: DS9 version 4.1\n"); + fileline.append("global color=green dashlist=8 3 width=1"); +// fileline.append("font=\"helvetica 10 normal\" select=1"); +// "highlite=1 dash=0 fixed=0 edit=1 move=1"+ +// "delete=1 include=1 source=1\n"; + ds9file << fileline; + + // Loop over regions in container + for (int i = 0; i < size(); ++i) { + ds9file << m_regions[i]->write() << "\n"; + } + + // Close file + ds9file.close(); + + // Store filename + m_filename = filename; + + } + + // ... otherwise, if file could not be opened then throw an exception + else { + throw GException::file_open_error(G_SAVE, filename); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print regions + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing region container information. + * + * Prints all regions into a string. + ***************************************************************************/ +std::string GSkyRegions::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GSkyRegions ==="); + + // Append information + result.append("\n"+gammalib::parformat("Number of regions")); + result.append(gammalib::str(size())); + + // Append regions + for (int i = 0; i < size(); ++i) { + result.append("\n"+m_regions[i]->print(chatter)); + } + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/***********************************************************************//** + * @brief Tells if direction is contained in one of the regions + * + * @param[in] dir A sky direction + * @return True or False + * + * Tells if direction is contained in one of the regions + ***************************************************************************/ + bool GSkyRegions::contains(const GSkyDir& dir) const + { + // Initialise return value + bool dir_is_in = false; + + // Loop over regions + for (int i = 0; i < size(); ++i) { + dir_is_in=m_regions[i]->contains(dir); + if (dir_is_in) break; + } + + // Return result + return dir_is_in; + } + + +/***********************************************************************//** + * @brief Tells if region overlaps one of the regions + * + * @param[in] reg A sky region + * @return True or False + * + * Tells if region overlaps one of the regions + ***************************************************************************/ + bool GSkyRegions::overlaps(const GSkyRegion& reg) const + { + // Initialise return value + bool reg_is_in = false; + + // Loop over regions + for (int i = 0; i < size(); ++i) { + reg_is_in=m_regions[i]->overlaps(reg); + if (reg_is_in) break; + } + + // Return result + return reg_is_in; + } + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GSkyRegions::init_members(void) +{ + // Initialise members + m_filename.clear(); + m_regions.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] regions region container. + * + * Makes a copy of all class members. All regions are deep copied, and the + * linear pointer array for parameter access through the GOptimizerPars + * base class is set. + ***************************************************************************/ +void GSkyRegions::copy_members(const GSkyRegions& regions) +{ + // Copy members + m_filename = regions.m_filename; + + // Copy regions + m_regions.clear(); + for (int i = 0; i < regions.m_regions.size(); ++i) { + m_regions.push_back((regions.m_regions[i]->clone())); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + * + * Deallocates all regions. The method loops over the region container and + * deallocates the memory that has been allocated before. + ***************************************************************************/ +void GSkyRegions::free_members(void) +{ + // Free regions + for (int i = 0; i < m_regions.size(); ++i) { + delete m_regions[i]; + m_regions[i] = NULL; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Return region index by name + * + * @param[in] name region name. + * @return region index (-1 if region name has not been found) + * + * Returns region index based on the specified @p name. If no region with the + * specified @p name is found the method returns -1. + ***************************************************************************/ +int GSkyRegions::get_index(const std::string& name) const +{ + // Initialise index + int index = -1; + + // Search region with specified name + for (int i = 0; i < size(); ++i) { + if (m_regions[i]->name() == name) { + index = i; + break; + } + } + + // Return index + return index; +} diff --git a/src/sky/Makefile.am b/src/sky/Makefile.am index db3a1d9..27f816d 100644 --- a/src/sky/Makefile.am +++ b/src/sky/Makefile.am @@ -9,6 +9,9 @@ INCLUDES = -I$(top_srcdir)/include \ sources = GSkyDir.cpp \ GSkyPixel.cpp \ GSkymap.cpp \ + GSkyRegions.cpp \ + GSkyRegion.cpp \ + GSkyRegionCircle.cpp \ GWcs.cpp \ GWcsRegistry.cpp \ GWcslib.cpp \ diff --git a/src/support/GTools.cpp b/src/support/GTools.cpp index 5e52579..d08fb56 100644 --- a/src/support/GTools.cpp +++ b/src/support/GTools.cpp @@ -824,3 +824,25 @@ bool gammalib::file_exists(const std::string& filename) // Return result return result; } + +/***********************************************************************//** + * @brief Checks if a substring is in a string . + * + * @param[in] str string you want to search in. + * @param[in] substring string you are looking for in str. + * + * Checks if substring is contained in str + ***************************************************************************/ +bool gammalib::contains(std::string str, std::string substring) +{ + // Initialise result + bool result = false; + + // checks if substring is in str + if (str.find(substring) != std::string::npos){ + result = true; + } + + // Return result + return result; +} diff --git a/src/xspec/GArf.cpp b/src/xspec/GArf.cpp new file mode 100644 index 0000000..68c9e1a --- /dev/null +++ b/src/xspec/GArf.cpp @@ -0,0 +1,481 @@ +/*************************************************************************** + * GArf.cpp - XSPEC Auxiliary Response File class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GArf.cpp + * @brief XSPEC Auxiliary Response File class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GArf.hpp" +#include "GException.hpp" +#include "GTools.hpp" +#include "GFitsBinTable.hpp" +#include "GFitsTableShortCol.hpp" +#include "GFitsTableFloatCol.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_AT "GArf::at(int&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GArf::GArf(void) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief File constructor + * + * @param[in] filename File name. + ***************************************************************************/ +GArf::GArf(const std::string& filename) +{ + // Initialise members + init_members(); + + // Load ARF file + load(filename); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Energy boundary constructor + * + * @param[in] ebds Energy boundaries. + ***************************************************************************/ +GArf::GArf(const GEbounds& ebds) +{ + // Initialise members + init_members(); + + // Set energy boundaries + m_ebounds = ebds; + + // Initialize spectral response + m_specresp.assign(ebds.size(), 0.0); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] arf Auxiliary Response File. + ***************************************************************************/ +GArf::GArf(const GArf& arf) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(arf); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GArf::~GArf(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] arf Auxiliary Response File. + * @return Auxiliary Response File. + ***************************************************************************/ +GArf& GArf::operator=(const GArf& arf) +{ + // Execute only if object is not identical + if (this != &arf) { + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(arf); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear object. + * + * Reset object to a clean initial state. + ***************************************************************************/ +void GArf::clear(void) +{ + // Free memory and initialise members + free_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone object + * + * @return Pointer to Auxiliary Response File. + ***************************************************************************/ +GArf* GArf::clone(void) const +{ + // Clone client + return new GArf(*this); +} + + +/***********************************************************************//** + * @brief Return content of spectral bin + * + * @param[in] index Bin index [0,...,size()-1]. + * + * @exception GException::out_of_range + * Spectral bin index is out of range. + * + * Returns reference to content of spectral bin with specified @p index. + ***************************************************************************/ +double& GArf::at(const int& index) +{ + // Raise exception if index is out of range + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_AT, index, 0, size()-1); + } + + // Return reference + return (m_specresp[index]); +} + + +/***********************************************************************//** + * @brief Return content of spectral bin (const version) + * + * @param[in] index Bin index [0,...,size()-1]. + * + * @exception GException::out_of_range + * Spectral bin index is out of range. + * + * Returns reference to content of spectral bin with specified @p index. + ***************************************************************************/ +const double& GArf::at(const int& index) const +{ + // Raise exception if index is out of range + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_AT, index, 0, size()-1); + } + + // Return reference + return (m_specresp[index]); +} + + +/***********************************************************************//** + * @brief Load Auxiliary Response File + * + * @param[in] filename File name. + ***************************************************************************/ +void GArf::load(const std::string& filename) +{ + // Clear any existing models + clear(); + + // Open FITS file + GFits file(filename); + + // Get ARF table + GFitsTable* table = file.table("SPECRESP"); + + // Read ARF data + read(table); + + // Close FITS file + file.close(); + + // Store filename + m_filename = filename; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save Auxiliary Response File + * + * @param[in] filename File name. + * @param[in] clobber Overwrite existing file? (defaults to true) + ***************************************************************************/ +void GArf::save(const std::string& filename, const bool& clobber) const +{ + // Open FITS file + GFits fits; + + // Write ARF into file + write(fits); + + // Close FITS file + fits.saveto(filename, clobber); + + // Store filename + m_filename = filename; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Read Auxiliary Response File + * + * @param[in] hdu ARF FITS table. + ***************************************************************************/ +void GArf::read(const GFitsTable* hdu) +{ + // 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; + } + + // Extract number of energy bins + int num = energy_lo->length(); + + // 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); + + } // endfor: looped over energy bins + + } // endif: HDU was valid + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write Auxiliary Response File + * + * @param[in] fits FITS file. + ***************************************************************************/ +void GArf::write(GFits& fits) const +{ + // Set column length + int length = size(); + + // Continue only if there are bins + if (length > 0) { + + // Create new binary table + GFitsBinTable* hdu = new GFitsBinTable; + + // Allocate floating point vector columns + GFitsTableFloatCol energy_lo("ENERG_LO", length); + GFitsTableFloatCol energy_hi("ENERG_HI", length); + GFitsTableFloatCol specresp("SPECRESP", length); + + // Fill columns + for (int i = 0; i < length; ++i) { + energy_lo(i) = m_ebounds.emin(i).keV(); + energy_hi(i) = m_ebounds.emax(i).keV(); + specresp(i) = m_specresp[i]; + } + + // Set column units + energy_lo.unit("keV"); + energy_hi.unit("keV"); + specresp.unit("cm^2"); + + // Set table attributes + hdu->extname("SPECRESP"); + + // Append columns to table + hdu->append_column(energy_lo); + hdu->append_column(energy_hi); + hdu->append_column(specresp); + + // Append HDU to FITS file + fits.append(*hdu); + + // Free binary table + delete hdu; + + } // endif: there were data to write + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print Auxiliary Response File + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing Auxiliary Response File information. + ***************************************************************************/ +std::string GArf::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GArf ==="); + + // Append energy boundary information + result.append("\n"+gammalib::parformat("Number of bins")); + result.append(gammalib::str(m_ebounds.size())); + result.append("\n"+gammalib::parformat("Energy range")); + result.append(m_ebounds.emin().print()); + result.append(" - "); + result.append(m_ebounds.emax().print()); + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GArf::init_members(void) +{ + // Initialise members + m_filename.clear(); + m_ebounds.clear(); + m_specresp.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] arf Auxiliary Response File. + ***************************************************************************/ +void GArf::copy_members(const GArf& arf) +{ + // Copy members + m_filename = arf.m_filename; + m_ebounds = arf.m_ebounds; + m_specresp = arf.m_specresp; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GArf::free_members(void) +{ + // Return + return; +} diff --git a/src/xspec/GPha.cpp b/src/xspec/GPha.cpp new file mode 100644 index 0000000..4284b34 --- /dev/null +++ b/src/xspec/GPha.cpp @@ -0,0 +1,571 @@ +/*************************************************************************** + * GPha.cpp - XSPEC Pulse Height Analyzer class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GPha.cpp + * @brief XSPEC Pulse Height Analyzer class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GPha.hpp" +#include "GException.hpp" +#include "GTools.hpp" +#include "GFitsBinTable.hpp" +#include "GFitsTableShortCol.hpp" +#include "GFitsTableFloatCol.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_AT "GPha::at(int&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GPha::GPha(void) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief File constructor + * + * @param[in] filename File name. + ***************************************************************************/ +GPha::GPha(const std::string& filename) +{ + // Initialise members + init_members(); + + // Load PHA file + load(filename); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Energy boundary constructor + * + * @param[in] ebds Energy boundaries. + ***************************************************************************/ +GPha::GPha(const GEbounds& ebds) +{ + // Initialise members + init_members(); + + // Set energy boundaries + m_ebounds = ebds; + + // Initialize spectrum + m_counts.assign(ebds.size(), 0.0); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] pha Pulse Height Analyzer spectrum. + ***************************************************************************/ +GPha::GPha(const GPha& pha) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(pha); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GPha::~GPha(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] pha Pulse Height Analyzer spectrum. + * @return Pulse Height Analyzer spectrum. + ***************************************************************************/ +GPha& GPha::operator= (const GPha& pha) +{ + // Execute only if object is not identical + if (this != &pha) { + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(pha); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear object. + * + * Reset object to a clean initial state. + ***************************************************************************/ +void GPha::clear(void) +{ + // Free memory and initialise members + free_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone object + * @return Pulse Height Analyzer spectrum. + ***************************************************************************/ +GPha* GPha::clone(void) const +{ + // Clone client + return new GPha(*this); +} + + +/***********************************************************************//** + * @brief Return content of spectral bin + * + * @param[in] index Bin index [0,...,size()-1]. + * + * @exception GException::out_of_range + * Spectral bin index is out of range. + * + * Returns reference to content of spectral bin with specified @p index. + ***************************************************************************/ +double& GPha::at(const int& index) +{ + // Raise exception if index is out of range + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_AT, index, 0, size()-1); + } + + // Return reference + return (m_counts[index]); +} + + +/***********************************************************************//** + * @brief Return content of spectral bin (const version) + * + * @param[in] index Bin index [0,...,size()-1]. + * + * @exception GException::out_of_range + * Spectral bin index is out of range. + * + * Returns reference to content of spectral bin with specified @p index. + ***************************************************************************/ +const double& GPha::at(const int& index) const +{ + // Raise exception if index is out of range + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_AT, index, 0, size()-1); + } + + // Return reference + return (m_counts[index]); +} + + +/***********************************************************************//** + * @brief Return index for a given energy + * + * Returns the number of counts in the spectrum. + ***************************************************************************/ +/* +int GPha::index(const GEnergy& energy) const +{ + // Initialise counts + double counts = 0.0; + + // Compute content + for (int i = 0; i < m_counts.size(); ++i) { + counts += m_counts[i]; + } + + // Return counts + return counts; +} +*/ + + +/***********************************************************************//** + * @brief Number of counts in spectrum + * + * Returns the number of counts in the spectrum. + ***************************************************************************/ +double GPha::counts(void) const +{ + // Initialise counts + double counts = 0.0; + + // Compute content + for (int i = 0; i < m_counts.size(); ++i) { + counts += m_counts[i]; + } + + // Return counts + return counts; +} + + +/***********************************************************************//** + * @brief Fill spectrum with a value. + * + * @param[in] energy Energy. + * @param[in] value Fill value (defaults to 1.0). + * + * Fills the specified @p value at a given @p energy in the spectrum. + ***************************************************************************/ +void GPha::fill(const GEnergy& energy, const double& value) +{ + // Get index + int index = m_ebounds.index(energy); + + // If index is not valid, then check for underflow, overflow or outflow + if (index == -1) { + if (energy < m_ebounds.emin()) { + m_underflow += value; + } + else if (energy >= m_ebounds.emax()) { + m_overflow += value; + } + else { + m_outflow += value; + } + } + + // ... otherwise fill the histogram + else { + m_counts[index] += value; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Load Pulse Height Analyzer spectrum + * + * @param[in] filename File name. + ***************************************************************************/ +void GPha::load(const std::string& filename) +{ + // Clear any existing models + clear(); + + // Open FITS file + GFits file(filename); + + // Read energy boundaries + m_ebounds.load(filename); + + // Get PHA table + GFitsTable* table = file.table("SPECTRUM"); + + // Read PHA data + read(table); + + // Close FITS file + file.close(); + + // Store filename + m_filename = filename; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save Pulse Height Analyzer spectrum + * + * @param[in] filename File name. + * @param[in] clobber Overwrite existing file? (defaults to true) + ***************************************************************************/ +void GPha::save(const std::string& filename, const bool& clobber) const +{ + // Open FITS file + GFits fits; + + // Write PHA into file + write(fits); + + // Close FITS file + fits.saveto(filename, clobber); + + // Store filename + m_filename = filename; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Read Pulse Height Analyzer spectrum + * + * @param[in] hdu PHA FITS table. + ***************************************************************************/ +void GPha::read(const GFitsTable* hdu) +{ + // Continue only if HDU is valid + if (hdu != NULL) { + + // Set expected column length + int length = m_ebounds.size(); + + // Initialize spectrum + m_counts.assign(length, 0.0); + + // Get pointer to data column + const GFitsTableCol* col_data = &(*hdu)["COUNTS"]; + + // Check whether column length is okay + //TODO + + // Copy data + for (int i = 0; i < length; ++i) { + m_counts[i] = col_data->real(i); + } + + } // endif: HDU was valid + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write Pulse Height Analyzer spectrum + * + * @param[in] fits FITS file. + ***************************************************************************/ +void GPha::write(GFits& fits) const +{ + // Set column length + int length = size(); + + // Continue only if there are bins + if (length > 0) { + + // Create new binary table + GFitsBinTable* hdu = new GFitsBinTable; + + // Allocate floating point vector columns + GFitsTableShortCol col_chan("CHANNEL", length); + GFitsTableFloatCol col_data("COUNTS", length); + GFitsTableFloatCol col_stat("STAT_ERR", length); + GFitsTableFloatCol col_syst("SYS_ERR", length); + GFitsTableShortCol col_qual("QUALITY", length); + GFitsTableShortCol col_grpg("GROUPING", length); + GFitsTableFloatCol col_area("AREASCAL", length); + GFitsTableFloatCol col_back("BACKSCAL", length); + + // Fill columns + for (int i = 0; i < length; ++i) { + col_chan(i) = i+1; // Channels start at 1 + col_data(i) = float(m_counts[i]); + } + + // Set table attributes + hdu->extname("SPECTRUM"); + + // Append columns to table + hdu->append_column(col_chan); + hdu->append_column(col_data); + hdu->append_column(col_stat); + hdu->append_column(col_syst); + hdu->append_column(col_qual); + hdu->append_column(col_grpg); + hdu->append_column(col_area); + hdu->append_column(col_back); + + // Append HDU to FITS file + fits.append(*hdu); + + // Free binary table + delete hdu; + + // Append energy boundaries + m_ebounds.write(&fits); + + } // endif: there were data to write + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print Pulse Height Analyzer spectrum + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing Pulse Height Analyzer spectrum information. + ***************************************************************************/ +std::string GPha::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GPha ==="); + + // Append energy boundary information + result.append("\n"+gammalib::parformat("Number of bins")); + result.append(gammalib::str(m_ebounds.size())); + result.append("\n"+gammalib::parformat("Energy range")); + result.append(m_ebounds.emin().print()); + result.append(" - "); + result.append(m_ebounds.emax().print()); + + // Append content information + result.append("\n"+gammalib::parformat("Total number of counts")); + result.append(gammalib::str(counts())); + result.append("\n"+gammalib::parformat("Underflow counts")); + result.append(gammalib::str(m_underflow)); + result.append("\n"+gammalib::parformat("Overflow counts")); + result.append(gammalib::str(m_overflow)); + result.append("\n"+gammalib::parformat("Outflow counts")); + result.append(gammalib::str(m_outflow)); + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GPha::init_members(void) +{ + // Initialise members + m_filename.clear(); + m_ebounds.clear(); + m_counts.clear(); + m_underflow = 0.0; + m_overflow = 0.0; + m_outflow = 0.0; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] pha Pulse Height Analyzer spectrum. + ***************************************************************************/ +void GPha::copy_members(const GPha& pha) +{ + // Copy members + m_filename = pha.m_filename; + m_ebounds = pha.m_ebounds; + m_counts = pha.m_counts; + m_underflow = pha.m_underflow; + m_overflow = pha.m_overflow; + m_outflow = pha.m_outflow; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GPha::free_members(void) +{ + // Return + return; +} diff --git a/src/xspec/GRmf.cpp b/src/xspec/GRmf.cpp new file mode 100644 index 0000000..0bfb300 --- /dev/null +++ b/src/xspec/GRmf.cpp @@ -0,0 +1,593 @@ +/*************************************************************************** + * GRmf.cpp - XSPEC Redistribution Matrix File class * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GRmf.cpp + * @brief XSPEC Redistribution Matrix File class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GRmf.hpp" +#include "GException.hpp" +#include "GTools.hpp" +#include "GFitsBinTable.hpp" +#include "GFitsTableShortCol.hpp" +#include "GFitsTableFloatCol.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_AT "GRmf::at(int&, int&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GRmf::GRmf(void) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief File constructor + * + * @param[in] filename File name. + ***************************************************************************/ +GRmf::GRmf(const std::string& filename) +{ + // Initialise members + init_members(); + + // Load RMF file + load(filename); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Energy boundary constructor + * + * @param[in] etrue True energy boundaries. + * @param[in] emeasured Measured energy boundaries. + ***************************************************************************/ +GRmf::GRmf(const GEbounds& etrue, const GEbounds& emeasured) +{ + // Initialise members + init_members(); + + // Set energy boundaries + m_ebds_true = etrue; + m_ebds_measured = emeasured; + + // Initialize matrix + m_matrix = GMatrixSparse(m_ebds_true.size(), m_ebds_measured.size()); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] rmf Redistribution Matrix File. + ***************************************************************************/ +GRmf::GRmf(const GRmf& rmf) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(rmf); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GRmf::~GRmf(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] rmf Redistribution Matrix File. + * @return Redistribution Matrix File. + ***************************************************************************/ +GRmf& GRmf::operator=(const GRmf& rmf) +{ + // Execute only if object is not identical + if (this != &rmf) { + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(rmf); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear object. + * + * Reset object to a clean initial state. + ***************************************************************************/ +void GRmf::clear(void) +{ + // Free memory and initialise members + free_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone object + * + * @return Redistribution Matrix File. + ***************************************************************************/ +GRmf* GRmf::clone(void) const +{ + // Clone object + return new GRmf(*this); +} + + +/***********************************************************************//** + * @brief Return content of redistribution matrix bin + * + * @param[in] itrue True energy index [0,...,ntrue()-1]. + * @param[in] imeasured Measured energy index [0,...,nmeasured()-1]. + * + * @exception GException::out_of_range + * Bin index is out of range. + * + * Returns reference to content of redistribution matrix bin bin with true + * energy index @p itrue and measured energy index @p imeasured. + ***************************************************************************/ +double& GRmf::at(const int& itrue, const int& imeasured) +{ + // Raise exception if indices are out of range + if (itrue < 0 || itrue >= ntrue()) { + throw GException::out_of_range(G_AT, itrue, 0, ntrue()-1); + } + // Raise exception if indices are out of range + if (imeasured < 0 || imeasured >= nmeasured()) { + throw GException::out_of_range(G_AT, imeasured, 0, nmeasured()-1); + } + + // Return reference + return (m_matrix(itrue, imeasured)); +} + + +/***********************************************************************//** + * @brief Return content of redistribution matrix bin (const version) + * + * @param[in] itrue True energy index [0,...,ntrue()-1]. + * @param[in] imeasured Measured energy index [0,...,nmeasured()-1]. + * + * @exception GException::out_of_range + * Bin index is out of range. + * + * Returns reference to content of redistribution matrix bin bin with true + * energy index @p itrue and measured energy index @p imeasured. + ***************************************************************************/ +const double& GRmf::at(const int& itrue, const int& imeasured) const +{ + // Raise exception if indices are out of range + if (itrue < 0 || itrue >= ntrue()) { + throw GException::out_of_range(G_AT, itrue, 0, ntrue()-1); + } + // Raise exception if indices are out of range + if (imeasured < 0 || imeasured >= nmeasured()) { + throw GException::out_of_range(G_AT, imeasured, 0, nmeasured()-1); + } + + // Return reference + return (m_matrix(itrue, imeasured)); +} + + +/***********************************************************************//** + * @brief Load Redistribution Matrix File + * + * @param[in] filename File name. + ***************************************************************************/ +void GRmf::load(const std::string& filename) +{ + // Clear any existing models + clear(); + + // Open FITS file + GFits file(filename); + + // Read measured energy boundaries + m_ebds_measured.load(filename); + + // Get RMF table + GFitsTable* table = file.table("MATRIX"); + + // Read RMF data + read(table); + + // Close FITS file + file.close(); + + // Store filename + m_filename = filename; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save Redistribution Matrix File + * + * @param[in] filename File name. + * @param[in] clobber Overwrite existing file? (defaults to true) + ***************************************************************************/ +void GRmf::save(const std::string& filename, const bool& clobber) const +{ + // Open FITS file + GFits fits; + + // Write RMF into file + write(fits); + + // Close FITS file + fits.saveto(filename, clobber); + + // Store filename + m_filename = filename; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Read Redistribution Matrix File + * + * @param[in] hdu RMF FITS table. + ***************************************************************************/ +void GRmf::read(const GFitsTable* hdu) +{ + // Continue only if HDU is valid + if (hdu != NULL) { + + // 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"]; + + // Set matrix rows and columns + int rows = energy_lo->length(); + int columns = m_ebds_measured.size(); + + // Initialize matrix + m_matrix = GMatrixSparse(rows, columns); + + // 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 + + } // endfor: looped over true energy bins + + } // endif: HDU was valid + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write Redistribution Matrix File + * + * @param[in] fits FITS file. + ***************************************************************************/ +void GRmf::write(GFits& fits) const +{ + // Set table length + int length = ntrue(); + + // Continue only if there are true energies + if (length > 0) { + + // Compute maximum number of matrix elements + int maxcolumns = 1; + for (int i = 0; i < length; i++) { + + // Search first non-zero element + int ifirst = 0; + for (; ifirst < nmeasured(); ++ifirst) { + if (m_matrix(i, ifirst) != 0.0) { + break; + } + } + + // Search last non-zero element + int ilast = nmeasured() - 1; + for (; ilast >= 0; --ilast) { + if (m_matrix(i, ilast) != 0.0) { + break; + } + } + + // Compute number of non-zero matrix elements + int num = ilast - ifirst + 1; + + // Set maximum + if (num > maxcolumns) { + maxcolumns = num; + } + + } // endfor: looped over true energies + + // Create new binary table + GFitsBinTable* hdu = new GFitsBinTable; + + // Allocate columns + GFitsTableFloatCol energy_lo("ENERG_LO", length); + GFitsTableFloatCol energy_hi("ENERG_HI", length); + GFitsTableShortCol n_grp("N_GRP", length); + GFitsTableShortCol f_chan("F_CHAN", length); + GFitsTableShortCol n_chan("N_CHAN", length); + GFitsTableFloatCol matrix("MATRIX", length, maxcolumns); + + // Fill columns + for (int itrue = 0; itrue < length; ++itrue) { + + // Search first non-zero element + int ifirst = 0; + for (; ifirst < nmeasured(); ++ifirst) { + if (m_matrix(itrue, ifirst) != 0.0) { + break; + } + } + + // Search last non-zero element + int ilast = nmeasured() - 1; + for (; ilast >= 0; --ilast) { + if (m_matrix(itrue, ilast) != 0.0) { + break; + } + } + + // Compute number of non-zero matrix elements + int num = ilast - ifirst + 1; + if (num < 0) { + ifirst = 0; + num = 0; + } + + // Set energy bin + energy_lo(itrue) = m_ebds_true.emin(itrue).keV(); + energy_hi(itrue) = m_ebds_true.emax(itrue).keV(); + + // Set group information + n_grp(itrue) = 1; + f_chan(itrue) = ifirst; + n_chan(itrue) = num; + + // Set matrix elements + for (int i = 0; i < num; ++i, ++ifirst) { + matrix(itrue, i) = m_matrix(itrue, ifirst); + } + + } // endfor: looped over true energy bins + + // Set column units + energy_lo.unit("keV"); + energy_hi.unit("keV"); + + // Set table attributes + hdu->extname("MATRIX"); + + // Append columns to table + hdu->append_column(energy_lo); + hdu->append_column(energy_hi); + hdu->append_column(n_grp); + hdu->append_column(f_chan); + hdu->append_column(n_chan); + hdu->append_column(matrix); + + // Append HDU to FITS file + fits.append(*hdu); + + // Free binary table + delete hdu; + + // Append measured energy boundaries + m_ebds_measured.write(&fits); + + } // endif: there were data to write + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print Redistribution Matrix File + * + * @param[in] chatter Chattiness (defaults to NORMAL). + * @return String containing Redistribution Matrix File information. + ***************************************************************************/ +std::string GRmf::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GRmf ==="); + + // Append energy boundary information + result.append("\n"+gammalib::parformat("Number of true energy bins")); + result.append(gammalib::str(m_ebds_true.size())); + result.append("\n"+gammalib::parformat("Number of measured bins")); + result.append(gammalib::str(m_ebds_measured.size())); + result.append("\n"+gammalib::parformat("True energy range")); + result.append(m_ebds_true.emin().print()); + result.append(" - "); + result.append(m_ebds_true.emax().print()); + result.append("\n"+gammalib::parformat("Measured energy range")); + result.append(m_ebds_measured.emin().print()); + result.append(" - "); + result.append(m_ebds_measured.emax().print()); + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GRmf::init_members(void) +{ + // Initialise members + m_filename.clear(); + m_ebds_true.clear(); + m_ebds_measured.clear(); + m_matrix.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] rmf Redistribution Matrix File. + ***************************************************************************/ +void GRmf::copy_members(const GRmf& rmf) +{ + // Copy members + m_filename = rmf.m_filename; + m_ebds_true = rmf.m_ebds_true; + m_ebds_measured = rmf.m_ebds_measured; + m_matrix = rmf.m_matrix; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GRmf::free_members(void) +{ + // Return + return; +} diff --git a/src/xspec/Makefile.am b/src/xspec/Makefile.am new file mode 100644 index 0000000..996a0f8 --- /dev/null +++ b/src/xspec/Makefile.am @@ -0,0 +1,16 @@ +## Process this file with automake to produce Makefile.in + +# Set include directories +INCLUDES = -I$(top_srcdir)/include \ + -I$(top_builddir)/src \ + -I$(top_srcdir)/src/support + +# Define sources for this directory +sources = GArf.cpp \ + GPha.cpp \ + GRmf.cpp + +# Build libtool library +noinst_LTLIBRARIES = libxspec.la +libxspec_la_SOURCES = $(sources) +libxspec_la_LDFLAGS = -version-info $(GAMMALIB_LT_VERSION) diff --git a/test/Makefile.am b/test/Makefile.am index a9ff828..530a952 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -71,6 +71,7 @@ check_PROGRAMS = test_GSupport \ test_GFits \ test_GXml \ test_GVO \ + test_GXspec \ test_GApplication \ test_GModel \ test_GSky \ @@ -92,6 +93,7 @@ TESTS = test_GSupport \ test_GFits \ test_GXml \ test_GVO \ + test_GXspec \ test_GApplication \ test_GModel \ test_GSky \ @@ -146,6 +148,11 @@ test_GVO_LDFLAGS = @LDFLAGS@ test_GVO_CPPFLAGS = @CPPFLAGS@ test_GVO_LDADD = $(top_srcdir)/src/libgamma.la # +test_GXspec_SOURCES = test_GXspec.cpp +test_GXspec_LDFLAGS = @LDFLAGS@ +test_GXspec_CPPFLAGS = @CPPFLAGS@ +test_GXspec_LDADD = $(top_srcdir)/src/libgamma.la +# test_GApplication_SOURCES = test_GApplication.cpp test_GApplication_LDFLAGS = @LDFLAGS@ test_GApplication_CPPFLAGS = @CPPFLAGS@ diff --git a/test/data/test_circle_radeczerozero.reg b/test/data/test_circle_radeczerozero.reg new file mode 100644 index 0000000..98fcb8e --- /dev/null +++ b/test/data/test_circle_radeczerozero.reg @@ -0,0 +1,4 @@ +# Region file format: DS9 version 4.1 +global color=green dashlist=8 3 width=1 font="helvetica 10 normal roman" select=1 highlite=1 dash=0 fixed=0 edit=1 move=1 delete=1 include=1 source=1 +fk5 +circle(0.,0.,10.) diff --git a/test/test_GSky.cpp b/test/test_GSky.cpp index 77cd96b..e3e6ccb 100644 --- a/test/test_GSky.cpp +++ b/test/test_GSky.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_GSky.cpp - test GSky classes * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2013 by Jurgen Knodlseder * + * copyright (C) 2010-2013 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -21,7 +21,7 @@ /** * @file test_GSky.cpp * @brief Testing of sky module. - * @author J. Knodlseder + * @author Juergen Knoedlseder */ /* __ Includes ___________________________________________________________ */ @@ -41,9 +41,11 @@ /* __ Debug definitions __________________________________________________ */ #define G_WCS_FORTH_BACK_PIXEL_DEBUG #define G_WCS_COPY_DEBUG + + /***********************************************************************//** -* @brief Set parameters and tests -***************************************************************************/ + * @brief Set parameters and tests + ***************************************************************************/ void TestGSky::set(void){ // Test name @@ -55,7 +57,9 @@ void TestGSky::set(void){ add_test(static_cast(&TestGSky::test_GSkymap_healpix_io),"Test Healpix GSkymap I/O"); add_test(static_cast(&TestGSky::test_GSkymap_wcs_construct),"Test WCS GSkymap constructors"); add_test(static_cast(&TestGSky::test_GSkymap_wcs_io),"Test WCS GSkymap I/O"); - + add_test(static_cast(&TestGSky::test_GSkyRegions_io),"Test the GSkyregion class - io"); + add_test(static_cast(&TestGSky::test_GSkyRegionCircle_construct),"Test the GSkyregionCircle class - constuct"); + add_test(static_cast(&TestGSky::test_GSkyRegionCircle_logic),"Test the GSkyregionCircle class - logic "); return; } @@ -353,7 +357,7 @@ void TestGSky::test_GWcslib(void) /*************************************************************************** - * Test: GSkymap_healpix_construct * + * @brief GSkymap_healpix_construct ***************************************************************************/ void TestGSky::test_GSkymap_healpix_construct(void) { @@ -499,7 +503,7 @@ void TestGSky::test_GSkymap_healpix_construct(void) /*************************************************************************** - * Test: GSkymap_healpix_io * + * @brief GSkymap_healpix_io ***************************************************************************/ void TestGSky::test_GSkymap_healpix_io(void) { @@ -567,7 +571,7 @@ void TestGSky::test_GSkymap_healpix_io(void) /*************************************************************************** - * Test: GSkymap_wcs_construct * + * @brief GSkymap_wcs_construct ***************************************************************************/ void TestGSky::test_GSkymap_wcs_construct(void) { @@ -626,7 +630,7 @@ void TestGSky::test_GSkymap_wcs_construct(void) /*************************************************************************** - * Test: GSkymap_wcs_io * + * @brief GSkymap_wcs_io ***************************************************************************/ void TestGSky::test_GSkymap_wcs_io(void) { @@ -700,7 +704,190 @@ void TestGSky::test_GSkymap_wcs_io(void) /*************************************************************************** - * Main test function * + * @brief GSkyRegionCircle_construct + ***************************************************************************/ +void TestGSky::test_GSkyRegionCircle_construct(void) +{ + // Define region for comparison + GSkyDir refdir_radeczerozero = GSkyDir(); + refdir_radeczerozero.radec_deg(0,0); + double refradius = 10.; + + // Test constructing: + test_try("Test constructor"); + try { + GSkyRegionCircle circle(refdir_radeczerozero,refradius); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + //Test constructing with radius 0 + test_try("Test constructor2"); + try { + GSkyRegionCircle circle(refdir_radeczerozero,-1); + test_try_failure(); + } + catch (GException::invalid_argument &e) { + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + //Test constructing with radius 0 + test_try("Test radius assingment after"); + try { + GSkyRegionCircle circle(refdir_radeczerozero,refradius); + circle.radius(-1.); + test_try_failure(); + } + catch (GException::invalid_argument &e) { + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Check radius assignment + GSkyRegionCircle refregion(refdir_radeczerozero,refradius); + double refradius_check = refregion.radius(); + test_value(refradius,refradius_check,1.0e-10, "Test radius assignment"); + + // Check solid angle assignment + double solidangle_check = refregion.solidangle(); + double solidangle = 2*gammalib::pi*(1- std::cos(refradius /180 * gammalib::pi)); + test_value(solidangle_check,solidangle,1.0e-10, "Test solid angle assignment"); + + //exit test + return; +} + +/*************************************************************************** + * @brief GSkyRegions_logic + ***************************************************************************/ +void TestGSky::test_GSkyRegionCircle_logic(void) +{ + + // set input files /filenames +// const std::string file1 = "test_GSkyRegionCircle1.reg"; +// const std::string file2 = "test_GSkyRegionCircle2.reg"; + + //set reference values + + GSkyDir refdir_radeczerozero = GSkyDir(); + refdir_radeczerozero.radec_deg(0,0); + + GSkyDir refdir_lbzerozero = GSkyDir(); + refdir_lbzerozero.lb_deg(0,0); + + GSkyDir refdir_outside_refregion = GSkyDir(); + refdir_outside_refregion.radec_deg(100,80); + + GSkyDir refdir_raoffset = GSkyDir(); + refdir_raoffset.radec_deg(5,0); + + GSkyDir refdir_decpole = GSkyDir(); + refdir_decpole.radec_deg(0,89); + + GSkyDir refdir_ndecpole = GSkyDir(); + refdir_ndecpole.radec_deg(100,89); + + GSkyDir refdir_rapole = GSkyDir(); + refdir_rapole.radec_deg(359,0); + + GSkyRegionCircle refregion(refdir_radeczerozero,10); + GSkyDir centre = refregion.centre(); + + GSkyRegionCircle refregion_smaller(refdir_radeczerozero,5); + GSkyRegionCircle refregion_larger(refdir_radeczerozero,20); + GSkyRegionCircle refregion_raoffset(refdir_raoffset,10); + GSkyRegionCircle refregion_rapole(refdir_rapole,3); + GSkyRegionCircle refregion_decpole(refdir_decpole,3); + + // Test contain dirs + test_assert(refregion.contains(refdir_radeczerozero),"test for containment"); + test_assert(refregion.contains(centre), "test if center in circle"); + test_assert(!refregion.contains(refdir_outside_refregion), "test2 for containment"); + + // Test contain regions + test_assert(refregion.contains(refregion_smaller),"test for containment region"); + test_assert(!refregion.contains(refregion_larger), "test for containment region2 "); + test_assert(refregion.contains(refregion), "test3 for containment region"); + + test_assert(refregion.contains(refdir_rapole), "rapole for containment region"); + test_assert(refregion_decpole.contains(refdir_ndecpole), "rapole for containment region"); + + // Test overlaps + test_assert(refregion.overlaps(refregion_smaller),"test for overlap"); + test_assert(refregion.overlaps(refregion_larger),"test2 for overlap"); + test_assert(refregion.overlaps(refregion_raoffset),"test3 for overlap"); + test_assert(!refregion.overlaps(refregion_decpole),"test4 for overlap"); + + // Exit test + return; +} + + +/*************************************************************************** + * @brief Test GSkyRegions input and output + ***************************************************************************/ +void TestGSky::test_GSkyRegions_io(void) +{ + // Set filenames + const std::string filereg_radeczerozero = "data/test_circle_radeczerozero.reg"; + + // Test regions loading + test_try("Test regions loading"); + try { + GSkyRegions regions; + regions.load(filereg_radeczerozero); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + +// //check if load does the proper thing +// GSkyRegions regs; +// GSkyRegion* reg; +// GSkyRegionCircle* circ; +// GSkyDir dir; +// double rad; +// +// regs.load(filereg_radeczerozero); +// reg = regs[0]; +// if (reg->type() == "circle"){ +// circ = (GSkyRegionCircle*)reg; +// dir = circ->centre(); +// rad = circ->radius(); +// } +// GSkyDir refdir_radeczerozero = GSkyDir(); +// refdir_radeczerozero.radec_deg(0,0); +// +// test_value(rad,10,1.0e-10, "Test radius assignment"); +// +// + + // Test regions saving + test_try("Test regions saving"); + try { + GSkyRegions regions; + regions.save("region.ds9"); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Exit test + return; +} + + +/*************************************************************************** + * @brief Main test function ***************************************************************************/ int main(void) { diff --git a/test/test_GSky.hpp b/test/test_GSky.hpp index 97e6c08..6a9b409 100644 --- a/test/test_GSky.hpp +++ b/test/test_GSky.hpp @@ -40,6 +40,9 @@ class TestGSky : public GTestSuite void test_GSkymap_healpix_io(void); void test_GSkymap_wcs_construct(void); void test_GSkymap_wcs_io(void); + void test_GSkyRegions_io(void); + void test_GSkyRegionCircle_construct(void); + void test_GSkyRegionCircle_logic(void); // Private methods private: diff --git a/test/test_GXspec.cpp b/test/test_GXspec.cpp new file mode 100644 index 0000000..f39df80 --- /dev/null +++ b/test/test_GXspec.cpp @@ -0,0 +1,369 @@ +/*************************************************************************** + * test_GXspec.hpp - Test Xspec module * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file test_GXspec.cpp + * @brief Implementation of unit tests for Xspec module + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "test_GXspec.hpp" +#include "GTools.hpp" + + +/***********************************************************************//** + * @brief Set parameters and tests + **************************************************************************/ +void TestGXspec::set(void) +{ + // Test name + name("GXspec"); + + // Append tests + append(static_cast(&TestGXspec::test_GPha), "Test GPha class"); + append(static_cast(&TestGXspec::test_GArf), "Test GArf class"); + append(static_cast(&TestGXspec::test_GRmf), "Test GRmf class"); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Test GPha class + **************************************************************************/ +void TestGXspec::test_GPha(void) +{ + // Test void constructor + test_try("GPha void constructor"); + try { + GPha pha; + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Test energy boundary constructor + test_try("GPha energy boundary constructor"); + try { + GEbounds ebds(10, GEnergy(0.1, "TeV"), GEnergy(10.0, "TeV")); + GPha pha(ebds); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Test filling and accessing + GEbounds ebds(9, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"), false); + GPha pha(ebds); + for (int i = 0; i < 12; i += 2) { + GEnergy energy(double(i), "TeV"); + pha.fill(energy); + } + test_value(pha.counts(), 4.0); + test_value(pha.underflow(), 1.0); + test_value(pha.overflow(), 1.0); + test_value(pha.outflow(), 0.0); + for (int i = 0; i < 9; i += 2) { + test_value(pha.at(i), 0.0); + test_value(pha[i], 0.0); + } + for (int i = 1; i < 9; i += 2) { + test_value(pha.at(i), 1.0); + test_value(pha[i], 1.0); + } + pha[0] = 5.0; + pha[1] = 3.7; + test_value(pha[0], 5.0); + test_value(pha[1], 3.7); + + // Test saving and loading + pha.save("pha.fits", true); + test_assert(pha.filename() == "pha.fits", + "Unexpected filename \""+pha.filename()+"\"."); + pha.load("pha.fits"); + test_assert(pha.filename() == "pha.fits", + "Unexpected filename \""+pha.filename()+"\"."); + test_value(pha[0], 5.0, 1.0e-6); + test_value(pha[1], 3.7, 1.0e-6); + test_value(pha.counts(), 11.7, 1.0e-6); + test_value(pha.underflow(), 0.0, 1.0e-6); + test_value(pha.overflow(), 0.0, 1.0e-6); + test_value(pha.outflow(), 0.0, 1.0e-6); + for (int i = 2; i < 9; i += 2) { + test_value(pha.at(i), 0.0); + test_value(pha[i], 0.0); + } + for (int i = 3; i < 9; i += 2) { + test_value(pha.at(i), 1.0); + test_value(pha[i], 1.0); + } + + // Test constructing + GPha pha2("pha.fits"); + test_assert(pha2.filename() == "pha.fits", + "Unexpected filename \""+pha2.filename()+"\"."); + test_value(pha2[0], 5.0, 1.0e-6); + test_value(pha2[1], 3.7, 1.0e-6); + test_value(pha2.counts(), 11.7, 1.0e-6); + test_value(pha2.underflow(), 0.0, 1.0e-6); + test_value(pha2.overflow(), 0.0, 1.0e-6); + test_value(pha2.outflow(), 0.0, 1.0e-6); + for (int i = 2; i < 9; i += 2) { + test_value(pha2.at(i), 0.0); + test_value(pha2[i], 0.0); + } + for (int i = 3; i < 9; i += 2) { + test_value(pha2.at(i), 1.0); + test_value(pha2[i], 1.0); + } + //std::cout << pha << std::endl; + //std::cout << pha2 << std::endl; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Test GArf class + **************************************************************************/ +void TestGXspec::test_GArf(void) +{ + // Test void constructor + test_try("GArf void constructor"); + try { + GArf arf; + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Test energy boundary constructor + test_try("GArf energy boundary constructor"); + try { + GEbounds ebds(10, GEnergy(0.1, "TeV"), GEnergy(10.0, "TeV")); + GArf arf(ebds); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Test filling and accessing + GEbounds ebds(9, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV")); + GArf arf(ebds); + for (int i = 0; i < 9; i += 2) { + arf[i] = 1.0; + } + for (int i = 0; i < 9; i += 2) { + test_value(arf.at(i), 1.0); + test_value(arf[i], 1.0); + } + for (int i = 1; i < 9; i += 2) { + test_value(arf.at(i), 0.0); + test_value(arf[i], 0.0); + } + arf[0] = 5.0; + arf[1] = 3.7; + test_value(arf[0], 5.0); + test_value(arf[1], 3.7); + + // Test saving and loading + arf.save("arf.fits", true); + test_assert(arf.filename() == "arf.fits", + "Unexpected filename \""+arf.filename()+"\"."); + arf.load("arf.fits"); + test_assert(arf.filename() == "arf.fits", + "Unexpected filename \""+arf.filename()+"\"."); + test_value(arf[0], 5.0, 1.0e-6); + test_value(arf[1], 3.7, 1.0e-6); + for (int i = 2; i < 9; i += 2) { + test_value(arf.at(i), 1.0); + test_value(arf[i], 1.0); + } + for (int i = 3; i < 9; i += 2) { + test_value(arf.at(i), 0.0); + test_value(arf[i], 0.0); + } + + // Test constructing + GArf arf2("arf.fits"); + test_assert(arf2.filename() == "arf.fits", + "Unexpected filename \""+arf2.filename()+"\"."); + test_value(arf2[0], 5.0, 1.0e-6); + test_value(arf2[1], 3.7, 1.0e-6); + for (int i = 2; i < 9; i += 2) { + test_value(arf2.at(i), 1.0); + test_value(arf2[i], 1.0); + } + for (int i = 3; i < 9; i += 2) { + test_value(arf2.at(i), 0.0); + test_value(arf2[i], 0.0); + } + //std::cout << arf << std::endl; + //std::cout << arf2 << std::endl; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Test GRmf class + **************************************************************************/ +void TestGXspec::test_GRmf(void) +{ + // Test void constructor + test_try("GRmf void constructor"); + try { + GRmf rmf; + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Test energy boundary constructor + test_try("GRmf energy boundary constructor"); + try { + GEbounds ebds(10, GEnergy(0.1, "TeV"), GEnergy(10.0, "TeV")); + GRmf rmf(ebds, ebds); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Test filling and accessing + GEbounds ebds(9, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV")); + GRmf rmf(ebds, ebds); + for (int i = 0; i < 9; ++i) { + for (int k = i; k < i+3 && k < 9; ++k) { + rmf(i, k) = 1.0; + } + } + for (int i = 0; i < 9; ++i) { + int k = 0; + for (; k < i; ++k) { + test_value(rmf.at(i,k), 0.0); + test_value(rmf(i,k), 0.0); + } + for (; k < i+3 && k < 9; ++k) { + test_value(rmf.at(i,k), 1.0); + test_value(rmf(i,k), 1.0); + } + for (; k < 9; ++k) { + test_value(rmf.at(i,k), 0.0); + test_value(rmf(i,k), 0.0); + } + } + + // Test saving and loading + rmf.save("rmf.fits", true); + test_assert(rmf.filename() == "rmf.fits", + "Unexpected filename \""+rmf.filename()+"\"."); + rmf.load("rmf.fits"); + test_assert(rmf.filename() == "rmf.fits", + "Unexpected filename \""+rmf.filename()+"\"."); + /* + for (int i = 0; i < 9; ++i) { + for (int k = 0; k < 9; ++k) { + std::cout << rmf(i,k) << " "; + } + std::cout << std::endl; + } + */ + for (int i = 0; i < 9; ++i) { + int k = 0; + for (; k < i; ++k) { + test_value(rmf.at(i,k), 0.0); + test_value(rmf(i,k), 0.0); + } + for (; k < i+3 && k < 9; ++k) { + test_value(rmf.at(i,k), 1.0); + test_value(rmf(i,k), 1.0); + } + for (; k < 9; ++k) { + test_value(rmf.at(i,k), 0.0); + test_value(rmf(i,k), 0.0); + } + } + + // Test constructing + GRmf rmf2("rmf.fits"); + test_assert(rmf2.filename() == "rmf.fits", + "Unexpected filename \""+rmf2.filename()+"\"."); + for (int i = 0; i < 9; ++i) { + int k = 0; + for (; k < i; ++k) { + test_value(rmf.at(i,k), 0.0); + test_value(rmf(i,k), 0.0); + } + for (; k < i+3 && k < 9; ++k) { + test_value(rmf.at(i,k), 1.0); + test_value(rmf(i,k), 1.0); + } + for (; k < 9; ++k) { + test_value(rmf.at(i,k), 0.0); + test_value(rmf(i,k), 0.0); + } + } + //std::cout << rmf << std::endl; + //std::cout << rmf2 << std::endl; + + // Return + return; +} + + +/*************************************************************************** + * @brief Main entry point for test executable + ***************************************************************************/ +int main(void) +{ + // Allocate test suit container + GTestSuites testsuites("Xspec module"); + + // Initially assume that we pass all tests + bool success = true; + + // Create a test suite + TestGXspec test; + + // Append test suite to the container + testsuites.append(test); + + // Run the testsuites + success = testsuites.run(); + + // Save test report + testsuites.save("reports/GXspec.xml"); + + // Return success status + return (success ? 0 : 1); +} diff --git a/test/test_GXspec.hpp b/test/test_GXspec.hpp new file mode 100644 index 0000000..99d5c99 --- /dev/null +++ b/test/test_GXspec.hpp @@ -0,0 +1,56 @@ +/*************************************************************************** + * test_GXspec.hpp - Test Xspec module * + * ----------------------------------------------------------------------- * + * copyright (C) 2013 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file test_GXspec.hpp + * @brief Definition of unit tests for Xspec module + * @author Juergen Knoedlseder + */ + +#ifndef TEST_GXSPEC_HPP +#define TEST_GXSPEC_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" + + +/***********************************************************************//** + * @class TestGXspec + * + * @brief Test suite for Xspec module + ***************************************************************************/ +class TestGXspec : public GTestSuite { + +public: + // Constructors and destructors + TestGXspec(void) : GTestSuite() {} + virtual ~TestGXspec(void) {} + + // Methods + virtual void set(void); + void test_GPha(void); + void test_GArf(void); + void test_GRmf(void); + +private: + // Private members +}; + +#endif /* TEST_GXSPEC_HPP */ diff --git a/test/test_GXspec.py b/test/test_GXspec.py new file mode 100755 index 0000000..7896e06 --- /dev/null +++ b/test/test_GXspec.py @@ -0,0 +1,118 @@ +# ========================================================================== +# This module performs unit tests for the GammaLib xspec module. +# +# Copyright (C) 2013 Juergen Knoedlseder +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +from gammalib import * +from math import * +import os + + +# ==================================== # +# Test class for GammaLib xspec module # +# ==================================== # +class Test(GPythonTestSuite): + """ + Test class for GammaLib xspec module. + """ + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + GPythonTestSuite.__init__(self) + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions. + """ + # Set test name + self.name("xspec") + + # Append tests + self.append(self.test_pha, "Test GPha") + self.append(self.test_arf, "Test GArf") + self.append(self.test_rmf, "Test GRmf") + + # Return + return + + # Test GPha class + def test_pha(self): + """ + Test GPha class. + """ + # Allocate GPha object and check its size + nmeasured = 10 + emeasured = GEbounds(nmeasured, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV")) + pha = GPha(emeasured) + self.test_value(pha.size(), nmeasured) + + # Set and retrieve a value + value = 3.7 + pha[3] = value + self.test_value(pha[3], value) + + # Return + return + + # Test GArf class + def test_arf(self): + """ + Test GArf class. + """ + # Allocate GArf object and check its size + nmeasured = 10 + emeasured = GEbounds(nmeasured, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV")) + arf = GArf(emeasured) + self.test_value(arf.size(), nmeasured) + + # Set and retrieve a value + value = 3.7 + arf[3] = value + self.test_value(arf[3], value) + + # Return + return + + # Test GRmf class + def test_rmf(self): + """ + Test GRmf class. + """ + # Allocate GRmf object and check its size + ntrue = 20 + nmeasured = 10 + etrue = GEbounds(ntrue, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV")) + emeasured = GEbounds(nmeasured, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV")) + rmf = GRmf(etrue, emeasured) + self.test_value(rmf.size(), ntrue*nmeasured) + self.test_value(rmf.ntrue(), ntrue) + self.test_value(rmf.nmeasured(), nmeasured) + + # Set and retrieve a value + value = 3.7 + rmf[3,5] = value + self.test_value(rmf[3,5], value) + + # Return + return diff --git a/test/test_python.py b/test/test_python.py index 4921d70..8a86baa 100755 --- a/test/test_python.py +++ b/test/test_python.py @@ -32,6 +32,7 @@ import test_GSky import test_GSupport import test_GTest import test_GXml +import test_GXspec # Try importing MWL tests try: @@ -73,17 +74,18 @@ if __name__ == '__main__': suites = GTestSuites("Python interface testing") # Allocate test suites and append them to the container - suite1 = test_GApplication.Test() - suite2 = test_GFits.Test() - suite3 = test_GLinalg.Test() - suite4 = test_GModel.Test() - suite5 = test_GNumerics.Test() - suite6 = test_GObservation.Test() - suite7 = test_GOptimizer.Test() - suite8 = test_GSky.Test() - suite9 = test_GSupport.Test() + suite1 = test_GApplication.Test() + suite2 = test_GFits.Test() + suite3 = test_GLinalg.Test() + suite4 = test_GModel.Test() + suite5 = test_GNumerics.Test() + suite6 = test_GObservation.Test() + suite7 = test_GOptimizer.Test() + suite8 = test_GSky.Test() + suite9 = test_GSupport.Test() suite10 = test_GTest.Test() suite11 = test_GXml.Test() + suite12 = test_GXspec.Test() # Setup unit tests suite1.set() @@ -97,6 +99,7 @@ if __name__ == '__main__': suite9.set() suite10.set() suite11.set() + suite12.set() # Append tests to container suites.append(suite1) @@ -110,6 +113,7 @@ if __name__ == '__main__': suites.append(suite9) suites.append(suite10) suites.append(suite11) + suites.append(suite12) # Optionally handle MWL suite if has_mwl: