Samuel Flis / gammalib

Commit 7656e4397eda628caeafbf50a1c2ca8fc6642e1c

Authored by Jürgen Knödlseder
2 parents fbc0e4ea 630f5f50

Merge branch 'classical-analysis' into integration 

Conflicts:
	.gitignore
	ChangeLog
	NEWS
Inline Side-by-side
Showing 52 changed files with 7760 additions and 38 deletions
.gitignore
  View file @7656e43
... ... @@ -54,7 +54,7 @@ aclocal.m4
54 54 autom4te.cache
55 55 config.guess
56 56 config.h
57   -config.h.in
  57 +config.h.in*
58 58 config.status
59 59 config.sub
60 60 configure
... ... @@ -65,6 +65,7 @@ test/test_url.dat
65 65 test/*.xml
66 66 test/reports/*.xml
67 67 test/*.trs
  68 +test/region.ds9
68 69 examples/python/*.xml
69 70 examples/python/*.html
70 71  
... ... @@ -82,9 +83,11 @@ test/test_GSky
82 83 test/test_GSupport
83 84 test/test_GVector
84 85 test/test_GXml
  86 +test/test_GXspec
85 87 test/test_GVO
86 88 test/test_LAT
87 89 test/test_MWL
88 90 test/test_CTA
89 91 test/test_COM
90 92 examples/cpp/readmodel/readmodel
  93 +
... ...
1 1 Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
2 2 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr)
3 3 Christoph Deil (Christoph.Deil@mpi-hd.mpg.de)
4   -Michael Mayer (michael.mayer@desy.de)
5 4 \ No newline at end of file
  5 +Michael Mayer (michael.mayer@desy.de)
  6 +Anneli Schulz (anneli.schulz@desy.de)
  7 +Chia-Chun Lu (Chia-Chun.Lu@mpi-hd.mpg.de)
6 8 \ No newline at end of file
... ...
ChangeLog
  View file @7656e43
1   -2013-09-07 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2013-11-01 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 0.8.0 released
4 4 ========================
... ... @@ -12,6 +12,9 @@
12 12 Optimize computations
13 13 Introduce gammalib namespace for constants and functions
14 14 GModelSpatialDiffuseCube reads ENERGIES extension
  15 + Add sky region handling
  16 + Add on/off analysis support for CTA
  17 + Add xspec interface support
15 18  
16 19 2012-12-18 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
17 20  
... ...
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 00-08-00
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -September, 7th, 2013
  5 +1 November 2013
6 6  
7 7  
8 8 1. Introduction
... ... @@ -221,6 +221,12 @@ Add a sin() and cos() pre computation cache to speed up repeated
221 221 distance computations using dist() or dist_deg(). Coordinate access
222 222 methods are now declared inline to provide also faster data access.
223 223  
  224 +Classes have been added to support handling of sky regions. The
  225 +GSkyRegion class is the abstract base class of all sky regions. The
  226 +GSkyRegionCircle class implements a circle on the sky. The DS9 region
  227 +format is used for loading and saving region information into a file.
  228 +Regions a gathered in a regions container, which is implement by the
  229 +GSkyRegions class.
224 230  
225 231 13. Support module
226 232 ------------------
... ... @@ -251,29 +257,41 @@ with a SAMP Hub. The method supports SAMP discovery and registering and
251 257 unregistering.
252 258  
253 259  
254   -16. Test module
  260 +16. Xspec module
  261 +----------------
  262 +A module has been added to provide an interface to the Xspec analysis
  263 +package. Three classes have been implemented that support the Xspec
  264 +file formats: GPha, GArf and GRmf.
  265 +
  266 +
  267 +17. Test module
255 268 ---------------
256 269 None
257 270  
258 271  
259   -17. Multi-wavelength interface
  272 +18. Multi-wavelength interface
260 273 ------------------------------
261 274 None
262 275  
263 276  
264   -18. Fermi-LAT interface
  277 +19. Fermi-LAT interface
265 278 ----------------------
266 279 None
267 280  
268 281  
269   -19. CTA interface
  282 +20. CTA interface
270 283 -----------------
271 284 The computations for diffuse map templates have been optimized by
272 285 introducing a pre computation cache for IRF and Npred values. This
273 286 considerably reduces the computing time for diffuse spatial maps.
274 287  
  288 +Support for on/off analysis has been added. The GCTAOnOffObservation
  289 +class gathers all information that is required for an on/off analysis
  290 +for a given run. The GCTAOnOffObservations container allows to analyse
  291 +several runs in one shot.
  292 +
275 293  
276   -20. COMPTEL interface
  294 +21. COMPTEL interface
277 295 ---------------------
278 296 None
279 297  
... ...
configure.ac
  View file @7656e43
... ... @@ -860,7 +860,8 @@ else
860 860 pyext/gammalib/sky_wrap.cpp pyext/gammalib/sky.py
861 861 pyext/gammalib/support_wrap.cpp pyext/gammalib/support.py
862 862 pyext/gammalib/test_wrap.cpp pyext/gammalib/test.py
863   - pyext/gammalib/xml_wrap.cpp pyext/gammalib/xml.py],,
  863 + pyext/gammalib/xml_wrap.cpp pyext/gammalib/xml.py
  864 + pyext/gammalib/xspec_wrap.cpp pyext/gammalib/xspec.py],,
864 865 [has_wrappers="no"])
865 866 if test "x$with_mwl" = "xyes"; then
866 867 AC_CHECK_FILES([pyext/gammalib/mwl_wrap.cpp pyext/gammalib/mwl.py],,
... ... @@ -1239,6 +1240,7 @@ AC_CONFIG_FILES([Makefile
1239 1240 src/fits/Makefile
1240 1241 src/xml/Makefile
1241 1242 src/vo/Makefile
  1243 + src/xspec/Makefile
1242 1244 src/sky/Makefile
1243 1245 src/opt/Makefile
1244 1246 src/obs/Makefile
... ...
examples/cpp/readmodel/Makefile
  View file @7656e43
... ... @@ -4,9 +4,5 @@ LDFLAGS=-L/usr/local/gamma/lib -lgamma
4 4 DEPS=
5 5 OBJ=readmodel.cpp
6 6  
7   -%o: *.cpp $(DEPS)
8   - $(CXX) -c -o $@ $< $(CFLAGS)
9   -
10 7 readmodel: $(OBJ)
11 8 $(CXX) -o $@ $^ $(CFLAGS) $(LDFLAGS)
12   -
... ...
include/GArf.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GArf.hpp - XSPEC Auxiliary Response File class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GArf.hpp
  23 + * @brief XSPEC Auxiliary Response File class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GARF_HPP
  28 +#define GARF_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GBase.hpp"
  33 +#include "GEbounds.hpp"
  34 +#include "GFits.hpp"
  35 +#include "GFitsTable.hpp"
  36 +
  37 +
  38 +/***********************************************************************//**
  39 + * @class GArf
  40 + *
  41 + * @brief Auxiliary Response File class
  42 + ***************************************************************************/
  43 +class GArf : public GBase {
  44 +
  45 +public:
  46 + // Constructors and destructors
  47 + GArf(void);
  48 + explicit GArf(const std::string& filename);
  49 + explicit GArf(const GEbounds& ebds);
  50 + GArf(const GArf& arf);
  51 + virtual ~GArf(void);
  52 +
  53 + // Operators
  54 + GArf& operator=(const GArf& arf);
  55 + double& operator[](const int& index);
  56 + const double& operator[](const int& index) const;
  57 +
  58 + // Methods
  59 + void clear(void);
  60 + GArf* clone(void) const;
  61 + int size(void) const;
  62 + double& at(const int& index);
  63 + const double& at(const int& index) const;
  64 + const GEbounds& ebounds(void) const;
  65 + void load(const std::string& filename);
  66 + void save(const std::string& filename, const bool& clobber = false) const;
  67 + void read(const GFitsTable* hdu);
  68 + void write(GFits& fits) const;
  69 + const std::string& filename(void) const;
  70 + std::string print(const GChatter& chatter = NORMAL) const;
  71 +
  72 +protected:
  73 + // Protected methods
  74 + void init_members(void);
  75 + void copy_members(const GArf& pha);
  76 + void free_members(void);
  77 +
  78 + // Protected members
  79 + mutable std::string m_filename; //!< Filename of origin
  80 + GEbounds m_ebounds; //!< Energy boundaries
  81 + std::vector<double> m_specresp; //!< Spectral response
  82 +};
  83 +
  84 +
  85 +/***********************************************************************//**
  86 + * @brief Return content of spectral bin
  87 + *
  88 + * @param[in] index Bin index [0,...,size()-1].
  89 + *
  90 + * Returns reference to content of spectral bin with specified @p index.
  91 + ***************************************************************************/
  92 +inline
  93 +double& GArf::operator[](const int& index)
  94 +{
  95 + return (m_specresp[index]);
  96 +}
  97 +
  98 +
  99 +/***********************************************************************//**
  100 + * @brief Return content of spectral bin (const version)
  101 + *
  102 + * @param[in] index Bin index [0,...,size()-1].
  103 + *
  104 + * Returns reference to content of spectral bin with specified @p index.
  105 + ***************************************************************************/
  106 +inline
  107 +const double& GArf::operator[](const int& index) const
  108 +{
  109 + return (m_specresp[index]);
  110 +}
  111 +
  112 +
  113 +/***********************************************************************//**
  114 + * @brief Return number of spectral bins
  115 + *
  116 + * @return Number of spectral bins.
  117 + *
  118 + * Returns the number of spectral bins.
  119 + ***************************************************************************/
  120 +inline
  121 +int GArf::size(void) const
  122 +{
  123 + return (m_specresp.size());
  124 +}
  125 +
  126 +
  127 +/***********************************************************************//**
  128 + * @brief Return energy boundaries
  129 + *
  130 + * @return Energy boundaries for all spectral bins.
  131 + *
  132 + * Returns the energy boundaries for all spectral bins.
  133 + ***************************************************************************/
  134 +inline
  135 +const GEbounds& GArf::ebounds(void) const
  136 +{
  137 + return m_ebounds;
  138 +}
  139 +
  140 +
  141 +/***********************************************************************//**
  142 + * @brief Return file name
  143 + *
  144 + * @return File name from which the ARF information has been read or into
  145 + * which ARF information has been saved.
  146 + *
  147 + * Returns the file name from which the ARF information has been read or into
  148 + * which ARF information has been saved. The returned string will be empty if
  149 + * no load() or save() method has been called before.
  150 + ***************************************************************************/
  151 +inline
  152 +const std::string& GArf::filename(void) const
  153 +{
  154 + return (m_filename);
  155 +}
  156 +
  157 +#endif /* GARF_HPP */
... ...
include/GPha.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GPha.hpp - XSPEC Pulse Height Analyzer class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GPha.hpp
  23 + * @brief XSPEC Pulse Height Analyzer class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GPHA_HPP
  28 +#define GPHA_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GBase.hpp"
  33 +#include "GEnergy.hpp"
  34 +#include "GEbounds.hpp"
  35 +#include "GFits.hpp"
  36 +#include "GFitsTable.hpp"
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class GPha
  41 + *
  42 + * @brief Pulse Height Analyzer class
  43 + ***************************************************************************/
  44 +class GPha : public GBase {
  45 +
  46 +public:
  47 + // Constructors and destructors
  48 + GPha(void);
  49 + explicit GPha(const std::string& filename);
  50 + explicit GPha(const GEbounds& ebds);
  51 + GPha(const GPha& pha);
  52 + virtual ~GPha(void);
  53 +
  54 + // Operators
  55 + GPha& operator=(const GPha& pha);
  56 + double& operator[](const int& index);
  57 + const double& operator[](const int& index) const;
  58 +
  59 + // Methods
  60 + void clear(void);
  61 + GPha* clone(void) const;
  62 + int size(void) const;
  63 + double& at(const int& index);
  64 + const double& at(const int& index) const;
  65 + const GEbounds& ebounds(void) const;
  66 + double counts(void) const;
  67 + const double& underflow(void) const;
  68 + const double& overflow(void) const;
  69 + const double& outflow(void) const;
  70 + void fill(const GEnergy& energy, const double& value = 1.0);
  71 + void load(const std::string& filename);
  72 + void save(const std::string& filename, const bool& clobber = false) const;
  73 + void read(const GFitsTable* hdu);
  74 + void write(GFits& fits) const;
  75 + const std::string& filename(void) const;
  76 + std::string print(const GChatter& chatter = NORMAL) const;
  77 +
  78 +protected:
  79 + // Protected methods
  80 + void init_members(void);
  81 + void copy_members(const GPha& pha);
  82 + void free_members(void);
  83 +
  84 + // Protected members
  85 + mutable std::string m_filename; //!< Filename of origin
  86 + GEbounds m_ebounds; //!< Energy boundaries
  87 + std::vector<double> m_counts; //!< Counts data
  88 + double m_underflow; //!< Number of underflowing events
  89 + double m_overflow; //!< Number of overflowing events
  90 + double m_outflow; //!< Number of outflowing events
  91 +};
  92 +
  93 +
  94 +/***********************************************************************//**
  95 + * @brief Return content of spectral bin
  96 + *
  97 + * @param[in] index Bin index [0,...,size()-1].
  98 + *
  99 + * Returns reference to content of spectral bin with specified @p index.
  100 + ***************************************************************************/
  101 +inline
  102 +double& GPha::operator[](const int& index)
  103 +{
  104 + return (m_counts[index]);
  105 +}
  106 +
  107 +
  108 +/***********************************************************************//**
  109 + * @brief Return content of spectral bin (const version)
  110 + *
  111 + * @param[in] index Bin index [0,...,size()-1].
  112 + *
  113 + * Returns reference to content of spectral bin with specified @p index.
  114 + ***************************************************************************/
  115 +inline
  116 +const double& GPha::operator[](const int& index) const
  117 +{
  118 + return (m_counts[index]);
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Return number of bins in spectrum
  124 + *
  125 + * @return Number of bins in spectrum.
  126 + *
  127 + * Returns the number of bins in the spectrum.
  128 + ***************************************************************************/
  129 +inline
  130 +int GPha::size(void) const
  131 +{
  132 + return (m_counts.size());
  133 +}
  134 +
  135 +
  136 +/***********************************************************************//**
  137 + * @brief Return energy boundaries
  138 + *
  139 + * @return Energy boundaries for all spectral bins.
  140 + *
  141 + * Returns the energy boundaries for all spectral bins.
  142 + ***************************************************************************/
  143 +inline
  144 +const GEbounds& GPha::ebounds(void) const
  145 +{
  146 + return m_ebounds;
  147 +}
  148 +
  149 +
  150 +/***********************************************************************//**
  151 + * @brief Return number of underflow counts
  152 + *
  153 + * @return Number of counts with energies lower than the first energy bin.
  154 + *
  155 + * Returns the number of counts with energies lower than the first energy
  156 + * bin.
  157 + ***************************************************************************/
  158 +inline
  159 +const double& GPha::underflow(void) const
  160 +{
  161 + return m_underflow;
  162 +}
  163 +
  164 +
  165 +/***********************************************************************//**
  166 + * @brief Return number of overflow counts
  167 + *
  168 + * @return Number of counts with energies larger than the last energy bin.
  169 + *
  170 + * Returns the number of counts with energies larger than the last energy
  171 + * bin.
  172 + ***************************************************************************/
  173 +inline
  174 +const double& GPha::overflow(void) const
  175 +{
  176 + return m_overflow;
  177 +}
  178 +
  179 +
  180 +/***********************************************************************//**
  181 + * @brief Return number of outflow counts
  182 + *
  183 + * @return Number of counts with energies between energy bins.
  184 + *
  185 + * Returns the number of counts with energies between energy bins.
  186 + ***************************************************************************/
  187 +inline
  188 +const double& GPha::outflow(void) const
  189 +{
  190 + return m_outflow;
  191 +}
  192 +
  193 +
  194 +/***********************************************************************//**
  195 + * @brief Return file name
  196 + *
  197 + * @return File name from which the PHA information has been read or into
  198 + * which PHA information has been saved.
  199 + *
  200 + * Returns the file name from which the PHA information has been read or into
  201 + * which PHA information has been saved. The returned string will be empty if
  202 + * no load() or save() method has been called before.
  203 + ***************************************************************************/
  204 +inline
  205 +const std::string& GPha::filename(void) const
  206 +{
  207 + return (m_filename);
  208 +}
  209 +
  210 +#endif /* GPHA_HPP */
... ...
include/GRmf.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GRmf.hpp - XSPEC Redistribution Matrix File class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GRmf.hpp
  23 + * @brief XSPEC Redistribution Matrix File class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GRMF_HPP
  28 +#define GRMF_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GBase.hpp"
  33 +#include "GEbounds.hpp"
  34 +#include "GMatrixSparse.hpp"
  35 +#include "GFits.hpp"
  36 +#include "GFitsTable.hpp"
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class GRmf
  41 + *
  42 + * @brief Redistribution Matrix File class
  43 + *
  44 + * Matrix rows are true energy, columns are channels
  45 + ***************************************************************************/
  46 +class GRmf : public GBase {
  47 +
  48 +public:
  49 + // Constructors and destructors
  50 + GRmf(void);
  51 + explicit GRmf(const std::string& filename);
  52 + explicit GRmf(const GEbounds& etrue, const GEbounds& emeasured);
  53 + GRmf(const GRmf& rmf);
  54 + virtual ~GRmf(void);
  55 +
  56 + // Operators
  57 + GRmf& operator=(const GRmf& rmf);
  58 + double& operator()(const int& itrue, const int& imeasured);
  59 + const double& operator()(const int& itrue, const int& imeasured) const;
  60 +
  61 + // Methods
  62 + void clear(void);
  63 + GRmf* clone(void) const;
  64 + int size(void) const;
  65 + int ntrue(void) const;
  66 + int nmeasured(void) const;
  67 + double& at(const int& itrue, const int& imeasured);
  68 + const double& at(const int& itrue, const int& imeasured) const;
  69 + const GEbounds& etrue(void) const;
  70 + const GEbounds& emeasured(void) const;
  71 + void load(const std::string& filename);
  72 + void save(const std::string& filename, const bool& clobber = false) const;
  73 + void read(const GFitsTable* hdu);
  74 + void write(GFits& fits) const;
  75 + const std::string& filename(void) const;
  76 + std::string print(const GChatter& chatter = NORMAL) const;
  77 +
  78 +protected:
  79 + // Protected methods
  80 + void init_members(void);
  81 + void copy_members(const GRmf& rmf);
  82 + void free_members(void);
  83 +
  84 + // Protected members
  85 + mutable std::string m_filename; //!< Filename of origin
  86 + GEbounds m_ebds_true; //!< True energy boundaries
  87 + GEbounds m_ebds_measured; //!< Measured energy boundaries
  88 + GMatrixSparse m_matrix; //!< Sparse redistribution matrix
  89 +
  90 +};
  91 +
  92 +
  93 +/***********************************************************************//**
  94 + * @brief Return content of redistribution matrix bin
  95 + *
  96 + * @param[in] itrue True energy index [0,...,ntrue()-1].
  97 + * @param[in] imeasured Measured energy index [0,...,nmeasured()-1].
  98 + *
  99 + * Returns reference to content of redistribution matrix bin bin with true
  100 + * energy index @p itrue and measured energy index @p imeasured.
  101 + ***************************************************************************/
  102 +inline
  103 +double& GRmf::operator()(const int& itrue, const int& imeasured)
  104 +{
  105 + return (m_matrix(itrue, imeasured));
  106 +}
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Return content of redistribution matrix bin (const version)
  111 + *
  112 + * @param[in] itrue True energy index [0,...,ntrue()-1].
  113 + * @param[in] imeasured Measured energy index [0,...,nmeasured()-1].
  114 + *
  115 + * Returns reference to content of redistribution matrix bin bin with true
  116 + * energy index @p itrue and measured energy index @p imeasured.
  117 + ***************************************************************************/
  118 +inline
  119 +const double& GRmf::operator()(const int& itrue, const int& imeasured) const
  120 +{
  121 + return (m_matrix(itrue, imeasured));
  122 +}
  123 +
  124 +
  125 +/***********************************************************************//**
  126 + * @brief Return number of redistribution matrix bins
  127 + *
  128 + * @return Number of redistribution matrix bins.
  129 + *
  130 + * Returns the number of redistribution matrix bins.
  131 + ***************************************************************************/
  132 +inline
  133 +int GRmf::size(void) const
  134 +{
  135 + return (m_matrix.rows()*m_matrix.columns());
  136 +}
  137 +
  138 +
  139 +/***********************************************************************//**
  140 + * @brief Return number of true energy bins in redistribution matrix
  141 + *
  142 + * @return Number true energy bins in redistribution matrix.
  143 + *
  144 + * Returns the number of true energy bins in redistribution matrix.
  145 + ***************************************************************************/
  146 +inline
  147 +int GRmf::ntrue(void) const
  148 +{
  149 + return (m_matrix.rows());
  150 +}
  151 +
  152 +
  153 +/***********************************************************************//**
  154 + * @brief Return number of measured energy bins in redistribution matrix
  155 + *
  156 + * @return Number measured energy bins in redistribution matrix.
  157 + *
  158 + * Returns the number of measured energy bins in redistribution matrix.
  159 + ***************************************************************************/
  160 +inline
  161 +int GRmf::nmeasured(void) const
  162 +{
  163 + return (m_matrix.columns());
  164 +}
  165 +
  166 +
  167 +/***********************************************************************//**
  168 + * @brief Return true energy boundaries
  169 + *
  170 + * @return True energy boundaries for redistribution matrix.
  171 + *
  172 + * Returns the true energy boundaries for redistribution matrix.
  173 + ***************************************************************************/
  174 +inline
  175 +const GEbounds& GRmf::etrue(void) const
  176 +{
  177 + return m_ebds_true;
  178 +}
  179 +
  180 +
  181 +/***********************************************************************//**
  182 + * @brief Return measured energy boundaries
  183 + *
  184 + * @return Measured energy boundaries for redistribution matrix.
  185 + *
  186 + * Returns the measured energy boundaries for redistribution matrix.
  187 + ***************************************************************************/
  188 +inline
  189 +const GEbounds& GRmf::emeasured(void) const
  190 +{
  191 + return m_ebds_measured;
  192 +}
  193 +
  194 +
  195 +/***********************************************************************//**
  196 + * @brief Return file name
  197 + *
  198 + * @return File name from which the RMF information has been read or into
  199 + * which RMF information has been saved.
  200 + *
  201 + * Returns the file name from which the RMF information has been read or into
  202 + * which RMF information has been saved. The returned string will be empty if
  203 + * no load() or save() method has been called before.
  204 + ***************************************************************************/
  205 +inline
  206 +const std::string& GRmf::filename(void) const
  207 +{
  208 + return (m_filename);
  209 +}
  210 +
  211 +#endif /* GRMF_HPP */
... ...
include/GSkyRegion.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegion.hpp - Abstract virtual sky region base class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GSkyRegion.hpp
  23 + * @brief Abstract sky region base class interface definition
  24 + * @author Pierrick Martin
  25 + */
  26 +
  27 +#ifndef GSKYREGION_HPP
  28 +#define GSKYREGION_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GBase.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GSkyPixel;
  37 +class GSkyDir;
  38 +
  39 +/***********************************************************************//**
  40 + * @class GSkyRegion
  41 + *
  42 + * @brief Abstract interface for the sky region class
  43 + *
  44 + * This class provides an abstract interface for a sky region. The sky region
  45 + * is defined by an array of parameters names and values specific to the
  46 + * derived class where the region type or shape is implemented.
  47 + *
  48 + * Accessible information elements are:
  49 + * - Type of the region (circle, rectangle,...)
  50 + * - Solid angle subtended by the region
  51 + *
  52 + * The object can be initialised from a string in the DS9 region file format,
  53 + * and return its description as a string in the DS9 region file format. The
  54 + * input/ouput to/from files is handled in the container class GSkyRegions.
  55 + ***************************************************************************/
  56 +class GSkyRegion : public GBase {
  57 +
  58 +public:
  59 + // Constructors and destructors
  60 + GSkyRegion(void);
  61 + GSkyRegion(const GSkyRegion& region);
  62 + virtual ~GSkyRegion(void);
  63 +
  64 + // Operators
  65 + virtual GSkyRegion& operator=(const GSkyRegion& region);
  66 +
  67 + // Pure virtual methods
  68 + virtual void clear(void) = 0;
  69 + virtual GSkyRegion* clone(void) const = 0;
  70 + virtual void read(const std::string& regstring) = 0;
  71 + virtual std::string write(void) const = 0;
  72 + virtual bool contains(const GSkyDir& dir) const = 0;
  73 + virtual bool overlaps(const GSkyRegion& reg) const = 0;
  74 + virtual bool contains(const GSkyRegion& reg) const = 0;
  75 + virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
  76 +
  77 + // Implemented methods
  78 + std::string type(void) const;
  79 + std::string name(void) const;
  80 + double solidangle(void) const;
  81 +
  82 +protected:
  83 + // Protected methods
  84 + void init_members(void);
  85 + void copy_members(const GSkyRegion& region);
  86 + void free_members(void);
  87 +
  88 + // Protected members
  89 + std::string m_type; //!< Type of the region (circle, rectangle,...)
  90 + std::string m_name; //!< Name of the region
  91 + double m_solid; //!< Solid angle subtended by the region
  92 +};
  93 +
  94 +
  95 +/***********************************************************************//**
  96 + * @brief Return region name
  97 + *
  98 + * @return region name
  99 + *
  100 + * Returns the region name.
  101 + ***************************************************************************/
  102 +inline
  103 +std::string GSkyRegion::name(void) const
  104 +{
  105 + return (m_name);
  106 +}
  107 +
  108 +/***********************************************************************//**
  109 + * @brief Return region type
  110 + *
  111 + * @return region type
  112 + *
  113 + * Returns the region type.
  114 + ***************************************************************************/
  115 +inline
  116 +std::string GSkyRegion::type(void) const
  117 +{
  118 + return (m_type);
  119 +}
  120 +
  121 +/***********************************************************************//**
  122 + * @brief Return solid angle
  123 + *
  124 + * @return Solid angle
  125 + *
  126 + * Returns the solid angle subtended by the region (in sr).
  127 + ***************************************************************************/
  128 +inline
  129 +double GSkyRegion::solidangle(void) const
  130 +{
  131 + return (m_solid);
  132 +}
  133 +
  134 +#endif /* GSKYREGION_HPP */
... ...
include/GSkyRegionCircle.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegionCircle.hpp - circular sky region class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GSkyRegionCircle.hpp
  23 + * @brief Abstract sky region base class interface definition
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +#ifndef GSKYREGIONCIRCLE_HPP
  28 +#define GSKYREGIONCIRCLE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GMath.hpp"
  34 +#include "GSkyDir.hpp"
  35 +#include "GSkyRegion.hpp"
  36 +
  37 +
  38 +/***********************************************************************//**
  39 + * @class GSkyRegionCircle
  40 + *
  41 + * @brief Abstract interface for the circular sky region class
  42 + *
  43 + * This class provides an implementation for a circular sky region. The sky
  44 + * region is defined by an array of parameters, the meaning of which is
  45 + * specific to the derived class where the region type or shape is defined.
  46 + * The parameters are GModelPar objects for convenience.
  47 + *
  48 + * The class holds several properties such as solid angle subtended by the
  49 + * region and computed through internal method compute_solid().
  50 + *
  51 + * To be clarified:
  52 + * - Do we want a member relating the region to an observation run ?
  53 + * - Constructor and read/write using XML may not be needed if we use DS9
  54 + * region file format ?
  55 + * - Replace GModelPar by double for the parameters (GModelPar is overkill) ?
  56 + *
  57 + ***************************************************************************/
  58 +class GSkyRegionCircle : public GSkyRegion {
  59 +
  60 +public:
  61 + // Constructors and destructors
  62 + GSkyRegionCircle(void);
  63 + GSkyRegionCircle(const GSkyRegionCircle& circle);
  64 + explicit GSkyRegionCircle(const std::string &line);
  65 + GSkyRegionCircle(GSkyDir &centre, const double &radius);
  66 + // GSkyRegionCircle(const double &ra, const double &dec, const double &radius);
  67 + virtual ~GSkyRegionCircle(void);
  68 +
  69 + // Operators
  70 + GSkyRegionCircle& operator=(const GSkyRegionCircle& region);
  71 +
  72 + // Implemented methods
  73 + void clear(void);
  74 + GSkyRegionCircle* clone(void) const;
  75 + GSkyDir centre(void) const;
  76 + double radius(void) const;
  77 + void radius(const double& radius);
  78 + void centre(const GSkyDir& centre);
  79 + void centre(const double& ra,const double& dec);
  80 + void read(const std::string& line);
  81 + std::string write(void) const;
  82 + bool contains(const GSkyDir& dir) const;
  83 + bool contains(const GSkyRegion& reg) const;
  84 + bool overlaps(const GSkyRegion& reg) const;
  85 + std::string print(const GChatter& chatter = NORMAL) const;
  86 +
  87 +protected:
  88 + // Protected methods
  89 + void init_members(void);
  90 + void copy_members(const GSkyRegionCircle& region);
  91 + void free_members(void);
  92 + void compute_solid(void);
  93 +
  94 + // Protected members
  95 + GSkyDir m_centre; //!< Centre or reference point of the region
  96 + double m_radius; //!< Radius of circular the region [deg]
  97 +};
  98 +
  99 +
  100 +/***********************************************************************//**
  101 + * @brief Return region type
  102 + *
  103 + * @return region type
  104 + *
  105 + * Returns the region type.
  106 + ***************************************************************************/
  107 +inline
  108 +double GSkyRegionCircle::radius(void) const
  109 +{
  110 + return (m_radius);
  111 +}
  112 +
  113 +/***********************************************************************//**
  114 + * @brief Return region centre
  115 + *
  116 + * @return region centre
  117 + *
  118 + * Returns the region centre.
  119 + ***************************************************************************/
  120 +inline
  121 +GSkyDir GSkyRegionCircle::centre(void) const
  122 +{
  123 + return (m_centre);
  124 +}
  125 +
  126 +
  127 +
  128 +/***********************************************************************//**
  129 + * @brief Set centre of region
  130 + *
  131 + * @param[in] dir Center of region.
  132 + ***************************************************************************/
  133 +inline
  134 +void GSkyRegionCircle::centre(const GSkyDir& dir)
  135 +{
  136 + // Set centre
  137 + m_centre = dir;
  138 +
  139 + // Return
  140 + return;
  141 +}
  142 +
  143 +
  144 +/***********************************************************************//**
  145 + * @brief Set centre values
  146 + *
  147 + * @param[in] ra Right ascension value.
  148 + * @param[in] dec declination value.
  149 + ***************************************************************************/
  150 +inline
  151 +void GSkyRegionCircle::centre(const double& ra, const double& dec)
  152 +{
  153 + // Set centre values
  154 + m_centre.radec_deg(ra,dec);
  155 +
  156 + // Return
  157 + return;
  158 +}
  159 +
  160 +#endif /* GSKYREGIONCIRCLE_HPP */
  161 +
... ...
include/GSkyRegions.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegions.hpp - Sky region container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GSkyRegions.hpp
  23 + * @brief Sky regions container class definition
  24 + * @author Pierrick Martin
  25 + */
  26 +
  27 +#ifndef GSKYREGIONS_HPP
  28 +#define GSKYREGIONS_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include <vector>
  33 +#include "GBase.hpp"
  34 +#include "GSkyRegion.hpp"
  35 +#include "GException.hpp"
  36 +#include "GSkyDir.hpp"
  37 +
  38 +/***********************************************************************//**
  39 + * @class GSkyRegions
  40 + *
  41 + * @brief Sky region container class
  42 + *
  43 + * This container class collects sky regions objects derived from the
  44 + * GSkyRegion abstract class. These can be accessed from the access
  45 + * operator[] or through the at() method with index checking. Access to
  46 + * the region can be done from index or name, so the name of each region
  47 + * in the container has to be unique.
  48 + *
  49 + * The object can be initialised from a DS9 region file, and has load and
  50 + * save methods from/to a DS9 region file.
  51 + *
  52 + ***************************************************************************/
  53 +class GSkyRegions : public GBase {
  54 +
  55 +public:
  56 + // Constructors and destructors
  57 + GSkyRegions(void);
  58 + GSkyRegions(const GSkyRegions& regions);
  59 + explicit GSkyRegions(const std::string& filename);
  60 + virtual ~GSkyRegions(void);
  61 +
  62 + // Operators
  63 + GSkyRegions& operator=(const GSkyRegions& regions);
  64 + GSkyRegion* operator[](const int& index);
  65 + const GSkyRegion* operator[](const int& index) const;
  66 + GSkyRegion* operator[](const std::string& name);
  67 + const GSkyRegion* operator[](const std::string& name) const;
  68 +
  69 + // Methods
  70 + void clear(void);
  71 + GSkyRegions* clone(void) const;
  72 + GSkyRegion* at(const int& index);
  73 + const GSkyRegion* at(const int& index) const;
  74 + int size(void) const;
  75 + bool isempty(void) const;
  76 + GSkyRegion* set(const int& index, const GSkyRegion& region);
  77 + GSkyRegion* set(const std::string& name, const GSkyRegion& region);
  78 + GSkyRegion* append(const GSkyRegion& region);
  79 + GSkyRegion* insert(const int& index, const GSkyRegion& region);
  80 + GSkyRegion* insert(const std::string& name, const GSkyRegion& region);
  81 + void remove(const int& index);
  82 + void remove(const std::string& name);
  83 + void reserve(const int& num);
  84 + void extend(const GSkyRegions& regions);
  85 + bool contains(const std::string& name) const;
  86 + bool contains(const GSkyDir& dir) const;
  87 + bool overlaps(const GSkyRegion& reg) const;
  88 + void load(const std::string& filename);
  89 + void save(const std::string& filename) const;
  90 + const std::string& filename(void) const;
  91 + std::string print(const GChatter& chatter = NORMAL) const;
  92 +
  93 +protected:
  94 + // Protected methods
  95 + void init_members(void);
  96 + void copy_members(const GSkyRegions& regions);
  97 + void free_members(void);
  98 + int get_index(const std::string& name) const;
  99 +
  100 + // Protected members
  101 + mutable std::string m_filename; //!< Filename of origin
  102 + std::vector<GSkyRegion*> m_regions; //!< List of regions
  103 +};
  104 +
  105 +
  106 +/***********************************************************************//**
  107 + * @brief Return pointer to region
  108 + *
  109 + * @param[in] index region index [0,...,size()-1].
  110 + *
  111 + * Returns a pointer to the region with the specified @p index.
  112 + ***************************************************************************/
  113 +inline
  114 +GSkyRegion* GSkyRegions::operator[](const int& index)
  115 +{
  116 + return (m_regions[index]);
  117 +}
  118 +
  119 +
  120 +/***********************************************************************//**
  121 + * @brief Return pointer to region (const version)
  122 + *
  123 + * @param[in] index region index [0,...,size()-1].
  124 + *
  125 + * Returns a const pointer to the region with the specified @p index.
  126 + ***************************************************************************/
  127 +inline
  128 +const GSkyRegion* GSkyRegions::operator[](const int& index) const
  129 +{
  130 + return (m_regions[index]);
  131 +}
  132 +
  133 +
  134 +/***********************************************************************//**
  135 + * @brief Return number of regions in container
  136 + *
  137 + * @return Number of regions in container.
  138 + *
  139 + * Returns the number of regions in the region container.
  140 + ***************************************************************************/
  141 +inline
  142 +int GSkyRegions::size(void) const
  143 +{
  144 + return (m_regions.size());
  145 +}
  146 +
  147 +
  148 +/***********************************************************************//**
  149 + * @brief Signals if there are no regions in container
  150 + *
  151 + * @return True if container is empty, false otherwise.
  152 + *
  153 + * Signals if the region container does not contain any region.
  154 + ***************************************************************************/
  155 +inline
  156 +bool GSkyRegions::isempty(void) const
  157 +{
  158 + return (m_regions.empty());
  159 +}
  160 +
  161 +
  162 +/***********************************************************************//**
  163 + * @brief Reserves space for regions in container
  164 + *
  165 + * @param[in] num Number of regions
  166 + *
  167 + * Reserves space for @p num regions in the container.
  168 + ***************************************************************************/
  169 +inline
  170 +void GSkyRegions::reserve(const int& num)
  171 +{
  172 + m_regions.reserve(num);
  173 + return;
  174 +}
  175 +
  176 +
  177 +/***********************************************************************//**
  178 + * @brief Return file name
  179 + *
  180 + * @return File name from which the regions have been loaded or into which
  181 + * the regions have been saved.
  182 + *
  183 + * Returns the file name from which the regions have been loaded or into
  184 + * which the regions have been saved. The returned string will be empty if
  185 + * no load() or save() method has been called before.
  186 + ***************************************************************************/
  187 +inline
  188 +const std::string& GSkyRegions::filename(void) const
  189 +{
  190 + return (m_filename);
  191 +}
  192 +
  193 +#endif /* GSKYREGIONS_HPP */
... ...
include/GTools.hpp
  View file @7656e43
... ... @@ -87,6 +87,7 @@ namespace gammalib {
87 87 bool file_exists(const std::string& filename);
88 88 bool isinfinite(const double& x);
89 89 bool isnotanumber(const double& x);
  90 + bool contains(std::string str1, std::string str2);
90 91 }
91 92  
92 93 /***********************************************************************//**
... ...
include/GammaLib.hpp
  View file @7656e43
... ... @@ -107,6 +107,11 @@
107 107 /* __ VO module __________________________________________________________ */
108 108 #include "GVOClient.hpp"
109 109  
  110 +/* __ XSPEC module _______________________________________________________ */
  111 +#include "GArf.hpp"
  112 +#include "GPha.hpp"
  113 +#include "GRmf.hpp"
  114 +
110 115 /* __ Application module _________________________________________________ */
111 116 #include "GApplication.hpp"
112 117 #include "GLog.hpp"
... ... @@ -132,6 +137,9 @@
132 137 #include "GSkyDir.hpp"
133 138 #include "GSkyPixel.hpp"
134 139 #include "GSkymap.hpp"
  140 +#include "GSkyRegions.hpp"
  141 +#include "GSkyRegion.hpp"
  142 +#include "GSkyRegionCircle.hpp"
135 143 #include "GWcs.hpp"
136 144 #include "GWcsRegistry.hpp"
137 145 #include "GWcslib.hpp"
... ...
include/Makefile.am
  View file @7656e43
... ... @@ -61,6 +61,9 @@ pkginclude_HEADERS = GammaLib.hpp \
61 61 GXmlAttribute.hpp \
62 62 GXmlPI.hpp \
63 63 GVOClient.hpp \
  64 + GPha.hpp \
  65 + GArf.hpp \
  66 + GRmf.hpp \
64 67 GApplication.hpp \
65 68 GLog.hpp \
66 69 GPars.hpp \
... ... @@ -75,6 +78,9 @@ pkginclude_HEADERS = GammaLib.hpp \
75 78 GSkyDir.hpp \
76 79 GSkyPixel.hpp \
77 80 GSkymap.hpp \
  81 + GSkyRegion.hpp \
  82 + GSkyRegions.hpp \
  83 + GSkyRegionCircle.hpp \
78 84 GWcs.hpp \
79 85 GWcsRegistry.hpp \
80 86 GWcslib.hpp \
... ...
inst/cta/Makefile.am
  View file @7656e43
... ... @@ -11,6 +11,8 @@ INCLUDES = -I$(top_srcdir)/include \
11 11 sources = src/GCTASupport.cpp \
12 12 src/GCTAException.cpp \
13 13 src/GCTAObservation.cpp \
  14 + src/GCTAOnOffObservation.cpp \
  15 + src/GCTAOnOffObservations.cpp \
14 16 src/GCTAEventList.cpp \
15 17 src/GCTAEventAtom.cpp \
16 18 src/GCTAEventCube.cpp \
... ... @@ -41,6 +43,8 @@ sources = src/GCTASupport.cpp \
41 43 # Define headers to be installed
42 44 pkginclude_HEADERS = include/GCTAException.hpp \
43 45 include/GCTAObservation.hpp \
  46 + include/GCTAOnOffObservation.hpp \
  47 + include/GCTAOnOffObservations.hpp \
44 48 include/GCTAEventList.hpp \
45 49 include/GCTAEventAtom.hpp \
46 50 include/GCTAEventCube.hpp \
... ...
inst/cta/include/GCTAOnOffObservation.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GCTAOnOffObservation.hpp - CTA on-off observation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Chia-Chun Lu & Christoph Deil *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservation.hpp
  23 + * @brief CTA on-off observation class definition
  24 + * @author Chia-Chun Lu & Christoph Deil
  25 + */
  26 +
  27 +#ifndef GCTAONOFFOBSERVATION_HPP
  28 +#define GCTAONOFFOBSERVATION_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GBase.hpp"
  33 +#include "GPha.hpp"
  34 +#include "GArf.hpp"
  35 +#include "GRmf.hpp"
  36 +#include "GCTAEventList.hpp"
  37 +#include "GCTAEventAtom.hpp"
  38 +#include "GCTAObservation.hpp"
  39 +#include "GSkyRegions.hpp"
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class GCTAOnOffObservation
  44 + *
  45 + * @brief CTA on-off observation class
  46 + ***************************************************************************/
  47 +class GCTAOnOffObservation : public GBase {
  48 +
  49 +public:
  50 + // Constructors and destructors
  51 + GCTAOnOffObservation(void);
  52 + GCTAOnOffObservation(const GEbounds& ereco, const GSkyRegions& on,
  53 + const GSkyRegions& off);
  54 + GCTAOnOffObservation(const GCTAOnOffObservation& obs);
  55 + virtual ~GCTAOnOffObservation(void);
  56 +
  57 + // Operators
  58 + GCTAOnOffObservation& operator=(const GCTAOnOffObservation& obs);
  59 +
  60 + // Methods
  61 + void clear(void);
  62 + GCTAOnOffObservation* clone(void) const;
  63 + void name(const std::string& name) { m_name = name; }
  64 + void instrument(const std::string& instrument) { m_instrument = instrument; }
  65 + void id(const std::string& id) { m_id = id; }
  66 + void off_regions(const GSkyRegions& regions) {m_off_regions = regions;}
  67 + void on_regions(const GSkyRegions& regions) {m_on_regions = regions;}
  68 + const std::string& name(void) const { return m_name; }
  69 + const std::string& instrument(void) const { return m_instrument; }
  70 + const std::string& id(void) const { return m_id; }
  71 + const GPha& on_spec(void) const { return m_on_spec; }
  72 + const GPha& off_spec(void) const { return m_off_spec; }
  73 + const GArf& arf(void) { return m_arf; };
  74 + const GRmf& rmf(void) {return m_rmf;}
  75 + void fill(const GCTAObservation& obs);
  76 + void compute_response(const GCTAObservation& obs,
  77 + const GEbounds& etrue);
  78 + void read(const GXmlElement& xml);
  79 + void write(GXmlElement& xml) const;
  80 + std::string print(const GChatter& chatter = NORMAL) const;
  81 +
  82 +protected:
  83 + // Protected methods
  84 + void init_members(void);
  85 + void copy_members(const GCTAOnOffObservation& obs);
  86 + void free_members(void);
  87 + void compute_arf(const GCTAObservation& obs);
  88 + void compute_rmf(const GCTAObservation& obs, const GEbounds& etrue);
  89 +
  90 + // Protected data members
  91 + std::string m_name; //!< Name
  92 + std::string m_instrument; //!< Instrument name
  93 + std::string m_id; //!< Observation identifier
  94 + GPha m_on_spec;
  95 + GPha m_off_spec;
  96 + GArf m_arf;
  97 + GRmf m_rmf;
  98 + GSkyRegions m_on_regions;
  99 + GSkyRegions m_off_regions;
  100 +};
  101 +
  102 +#endif /* GCTAONOFFOBSERVATION_HPP */
... ...
inst/cta/include/GCTAOnOffObservations.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GCTAOnOffObservations.hpp - CTA on-off observation container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservations.hpp
  23 + * @brief CTA on-off observation container class definition
  24 + * @author Pierrick Martin
  25 + */
  26 +
  27 +#ifndef GCTAONOFFOBSERVATIONS_HPP
  28 +#define GCTAONOFFOBSERVATIONS_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include <vector>
  33 +#include "GContainer.hpp"
  34 +#include "GCTAOnOffObservation.hpp"
  35 +#include "GModels.hpp"
  36 +
  37 +
  38 +/***********************************************************************//**
  39 + * @class GCTAOnOffObservations
  40 + *
  41 + * @brief ON/OFF Observation container class
  42 + *
  43 + * This class in a container of GCTAOnOffObservation objects. Still some
  44 + * work to do and things to be clarified...
  45 + ***************************************************************************/
  46 +class GCTAOnOffObservations : public GContainer {
  47 +
  48 +public:
  49 + // Constructors and destructors
  50 + GCTAOnOffObservations(void);
  51 + GCTAOnOffObservations(const GCTAOnOffObservations& obs);
  52 + explicit GCTAOnOffObservations(const std::string& filename);
  53 + virtual ~GCTAOnOffObservations(void);
  54 +
  55 + // Operators
  56 + GCTAOnOffObservations& operator=(const GCTAOnOffObservations& obs);
  57 + GCTAOnOffObservation* operator[](const int& index);
  58 + const GCTAOnOffObservation* operator[](const int& index) const;
  59 +
  60 + // Methods
  61 + void clear(void);
  62 + GCTAOnOffObservations* clone(void) const;
  63 + GCTAOnOffObservation* at(const int& index);
  64 + const GCTAOnOffObservation* at(const int& index) const;
  65 + int size(void) const;
  66 + bool isempty(void) const;
  67 + GCTAOnOffObservation* set(const int& index, const GCTAOnOffObservation& obs);
  68 + GCTAOnOffObservation* append(const GCTAOnOffObservation& obs);
  69 + GCTAOnOffObservation* insert(const int& index, const GCTAOnOffObservation& obs);
  70 + void remove(const int& index);
  71 + void reserve(const int& num);
  72 + void extend(const GCTAOnOffObservations& obs);
  73 + bool contains(const std::string& instrument,
  74 + const std::string& id) const;
  75 + void load(const std::string& filename);
  76 + void save(const std::string& filename) const;
  77 + void read(const GXml& xml);
  78 + void write(GXml& xml) const;
  79 + void models(const GModels& models);
  80 + void models(const std::string& filename);
  81 + const GModels& models(void) const;
  82 + std::string print(const GChatter& chatter = NORMAL) const;
  83 +
  84 + /*
  85 + // To do
  86 + Implement fit functions for both cases summed over runs or not ?
  87 + Implement significance or have it at ctool level ?
  88 + Implement flux points computation and fit of a spectral function ?
  89 + // Do we want that ? If yes, need to retrieve inline functions from GObservations...
  90 + double npred(void) const;
  91 + // Left aside at the moment
  92 + void optimize(GOptimizer& opt);
  93 + */
  94 +
  95 +protected:
  96 + // Protected methods
  97 + void init_members(void);
  98 + void copy_members(const GCTAOnOffObservations& obs);
  99 + void free_members(void);
  100 + int get_index(const std::string& instrument,
  101 + const std::string& id) const;
  102 +
  103 + // Protected members
  104 + std::vector<GCTAOnOffObservation*> m_obs; //!< List of observations
  105 + GModels m_models; //!< List of models
  106 +};
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Return pointer to observation
  111 + *
  112 + * @param[in] index Observation index [0,...,size()-1].
  113 + * @return Observation.
  114 + *
  115 + * Returns a pointer to the observation with specified @p index.
  116 + ***************************************************************************/
  117 +inline
  118 +GCTAOnOffObservation* GCTAOnOffObservations::operator[](const int& index)
  119 +{
  120 + return m_obs[index];
  121 +}
  122 +
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Return pointer to observation (const version)
  126 + *
  127 + * @param[in] index Observation index [0,...,size()-1].
  128 + *
  129 + * Returns a const pointer to the observation with specified @p index.
  130 + ***************************************************************************/
  131 +inline
  132 +const GCTAOnOffObservation* GCTAOnOffObservations::operator[](const int& index) const
  133 +{
  134 + return m_obs[index];
  135 +}
  136 +
  137 +
  138 +/***********************************************************************//**
  139 + * @brief Return number of observations in container
  140 + *
  141 + * @return Number of observations in container.
  142 + *
  143 + * Returns the number of observations in the observation container.
  144 + ***************************************************************************/
  145 +inline
  146 +int GCTAOnOffObservations::size(void) const
  147 +{
  148 + return (m_obs.size());
  149 +}
  150 +
  151 +
  152 +/***********************************************************************//**
  153 + * @brief Signals if there are no observations in container
  154 + *
  155 + * @return True if container is empty, false otherwise.
  156 + *
  157 + * Signals if the observation container does not contain any observation.
  158 + ***************************************************************************/
  159 +inline
  160 +bool GCTAOnOffObservations::isempty(void) const
  161 +{
  162 + return (m_obs.empty());
  163 +}
  164 +
  165 +
  166 +/***********************************************************************//**
  167 + * @brief Reserves space for observations in container
  168 + *
  169 + * @param[in] num Number of observations.
  170 + *
  171 + * Reserves space for @p num observations in the container.
  172 + ***************************************************************************/
  173 +inline
  174 +void GCTAOnOffObservations::reserve(const int& num)
  175 +{
  176 + m_obs.reserve(num);
  177 + return;
  178 +}
  179 +
  180 +
  181 +/***********************************************************************//**
  182 + * @brief Set model container
  183 + *
  184 + * @param[in] models Model container.
  185 + *
  186 + * Sets the model container for the observations.
  187 + ***************************************************************************/
  188 +inline
  189 +void GCTAOnOffObservations::models(const GModels& models)
  190 +{
  191 + m_models = models;
  192 + return;
  193 +}
  194 +
  195 +
  196 +/***********************************************************************//**
  197 + * @brief Return model container
  198 + *
  199 + * @return Model container.
  200 + *
  201 + * Returns the model container of the observations.
  202 + ***************************************************************************/
  203 +inline
  204 +const GModels& GCTAOnOffObservations::models(void) const
  205 +{
  206 + return m_models;
  207 +}
  208 +
  209 +#endif /* GCTAONOFFOBSERVATIONS_HPP */
... ...
inst/cta/pyext/GCTAOnOffObservation.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GCTAOnOffObservation.i - CTA on-off observation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservation.i
  23 + * @brief CTA on-off observation class definition
  24 + * @author Michael Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCTAOnOffObservation.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCTAOnOffObservation
  34 + *
  35 + * @brief CTA on-off observation class
  36 + ***************************************************************************/
  37 +class GCTAOnOffObservation : public GBase {
  38 +public:
  39 + // Constructors and destructors
  40 + GCTAOnOffObservation(void);
  41 + GCTAOnOffObservation(const GEbounds& ereco, const GSkyRegions& on,
  42 + const GSkyRegions& off);
  43 + GCTAOnOffObservation(const GCTAOnOffObservation& obs);
  44 + virtual ~GCTAOnOffObservation(void);
  45 +
  46 + // Methods
  47 + void clear(void);
  48 + GCTAOnOffObservation* clone(void) const;
  49 + void name(const std::string& name) { m_name = name; }
  50 + void instrument(const std::string& instrument) { m_instrument = instrument; }
  51 + void id(const std::string& id) { m_id = id; }
  52 + void off_regions(const GSkyRegions& regions) {m_off_regions = regions;}
  53 + void on_regions(const GSkyRegions& regions) {m_on_regions = regions;}
  54 + const std::string& name(void) const { return m_name; }
  55 + const std::string& instrument(void) const { return m_instrument; }
  56 + const std::string& id(void) const { return m_id; }
  57 + const GPha& on_spec(void) const { return m_on_spec; }
  58 + const GPha& off_spec(void) const { return m_off_spec; }
  59 + const GArf& arf(void) { return m_arf; };
  60 + const GRmf& rmf(void) {return m_rmf;}
  61 + void fill(const GCTAObservation& obs);
  62 + void compute_response(const GCTAObservation& obs,
  63 + const GEbounds& etrue);
  64 + void read(const GXmlElement& xml);
  65 + void write(GXmlElement& xml) const;
  66 +};
  67 +
  68 +
  69 +/***********************************************************************//**
  70 + * @brief GCTOnOffAObservation class extension
  71 + ***************************************************************************/
  72 +%extend GCTAOnOffObservation {
  73 + GCTAOnOffObservation copy() {
  74 + return (*self);
  75 + }
  76 +}
... ...
inst/cta/pyext/GCTAOnOffObservations.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GCTAOnOffObservations.i - CTA on-off observation container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservations.i
  23 + * @brief CTA on-off observation container class definition
  24 + * @author Pierrick Martin
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCTAOnOffObservations.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCTAOnOffObservations
  34 + *
  35 + * @brief ON/OFF Observation container class
  36 + ***************************************************************************/
  37 +class GCTAOnOffObservations : public GContainer {
  38 +public:
  39 + // Constructors and destructors
  40 + GCTAOnOffObservations(void);
  41 + GCTAOnOffObservations(const GCTAOnOffObservations& obs);
  42 + explicit GCTAOnOffObservations(const std::string& filename);
  43 + virtual ~GCTAOnOffObservations(void);
  44 +
  45 + // Methods
  46 + void clear(void);
  47 + GCTAOnOffObservations* clone(void) const;
  48 + GCTAOnOffObservation* at(const int& index);
  49 + const GCTAOnOffObservation* at(const int& index) const;
  50 + int size(void) const;
  51 + bool isempty(void) const;
  52 + GCTAOnOffObservation* set(const int& index, const GCTAOnOffObservation& obs);
  53 + GCTAOnOffObservation* append(const GCTAOnOffObservation& obs);
  54 + GCTAOnOffObservation* insert(const int& index, const GCTAOnOffObservation& obs);
  55 + void remove(const int& index);
  56 + void reserve(const int& num);
  57 + void extend(const GCTAOnOffObservations& obs);
  58 + bool contains(const std::string& instrument,
  59 + const std::string& id) const;
  60 + void load(const std::string& filename);
  61 + void save(const std::string& filename) const;
  62 + void read(const GXml& xml);
  63 + void write(GXml& xml) const;
  64 + void models(const GModels& models);
  65 + void models(const std::string& filename);
  66 + const GModels& models(void) const;
  67 +};
  68 +
  69 +
  70 +/***********************************************************************//**
  71 + * @brief GCTAOnOffObservations class extension
  72 + ***************************************************************************/
  73 +%extend GCTAOnOffObservations {
  74 + GCTAOnOffObservation* __getitem__(const int& index) {
  75 + if (index >= 0 && index < self->size()) {
  76 + return (*self)[index];
  77 + }
  78 + else {
  79 + throw GException::out_of_range("__getitem__(int)", index, self->size());
  80 + }
  81 + }
  82 + void __setitem__(const int& index, const GCTAOnOffObservation& val) {
  83 + if (index>=0 && index < self->size()) {
  84 + self->set(index, val);
  85 + return;
  86 + }
  87 + else {
  88 + throw GException::out_of_range("__setitem__(int)", index, self->size());
  89 + }
  90 + }
  91 + int __len__() {
  92 + return (self->size());
  93 + }
  94 + GCTAOnOffObservations copy() {
  95 + return (*self);
  96 + }
  97 +};
... ...
inst/cta/pyext/cta.i
  View file @7656e43
... ... @@ -63,6 +63,8 @@
63 63  
64 64 /* __ CTA ________________________________________________________________ */
65 65 %include "GCTAObservation.i"
  66 +%include "GCTAOnOffObservation.i"
  67 +%include "GCTAOnOffObservations.i"
66 68 %include "GCTAEventCube.i"
67 69 %include "GCTAEventList.i"
68 70 %include "GCTAEventBin.i"
... ...
inst/cta/src/GCTAOnOffObservation.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GCTAOnOffObservation.cpp - CTA on-off observation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Chia-Chun Lu & Christoph Deil *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservation.cpp
  23 + * @brief CTA on-off observation class implementation
  24 + * @author Chia-Chun Lu & Christoph Deil
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GCTAOnOffObservation.hpp"
  32 +#include "GTools.hpp"
  33 +
  34 +/* __ Method name definitions ____________________________________________ */
  35 +#define G_WRITE "GCTAOnOffObservation::write(GXmlElement&)"
  36 +#define G_READ "GCTAOnOffObservation::read(GXmlElement&)"
  37 +#define G_FILL "GCTAOnOffObservation::fill(GCTAObservation&)"
  38 +
  39 +/* __ Macros _____________________________________________________________ */
  40 +
  41 +/* __ Coding definitions _________________________________________________ */
  42 +
  43 +/* __ Debug definitions __________________________________________________ */
  44 +
  45 +
  46 +/*==========================================================================
  47 + = =
  48 + = Constructors/destructors =
  49 + = =
  50 + ==========================================================================*/
  51 +
  52 +/***********************************************************************//**
  53 + * @brief Void constructor
  54 + ***************************************************************************/
  55 +GCTAOnOffObservation::GCTAOnOffObservation(void)
  56 +{
  57 + // Initialise private members
  58 + init_members();
  59 +
  60 + // Return
  61 + return;
  62 +}
  63 +
  64 +
  65 +/***********************************************************************//**
  66 + * @brief Copy constructor
  67 + *
  68 + * @param[in] obs CTA on-off observation.
  69 + ***************************************************************************/
  70 +GCTAOnOffObservation::GCTAOnOffObservation(const GCTAOnOffObservation& obs)
  71 +{
  72 + // Initialise private
  73 + init_members();
  74 +
  75 + // Copy members
  76 + copy_members(obs);
  77 +
  78 + // Return
  79 + return;
  80 +}
  81 +
  82 +
  83 +/***********************************************************************//**
  84 + * @brief Parameter constructor
  85 + *
  86 + * @param[in] ereco Reconstructed energy bins.
  87 + * @param[in] on On regions.
  88 + * @param[in] off Off regions.
  89 + ***************************************************************************/
  90 +GCTAOnOffObservation::GCTAOnOffObservation(const GEbounds& ereco,
  91 + const GSkyRegions& on,
  92 + const GSkyRegions& off)
  93 +{
  94 + // Initialise private
  95 + init_members();
  96 +
  97 + // Initialise spectra
  98 + m_on_spec = GPha(ereco);
  99 + m_off_spec = GPha(ereco);
  100 +
  101 + // Store regions
  102 + m_on_regions = on;
  103 + m_off_regions = off;
  104 +
  105 + // Return
  106 + return;
  107 +}
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Destructor
  112 + ***************************************************************************/
  113 +GCTAOnOffObservation::~GCTAOnOffObservation(void)
  114 +{
  115 + // Free members
  116 + free_members();
  117 +
  118 + // Return
  119 + return;
  120 +}
  121 +
  122 +
  123 +/*==========================================================================
  124 + = =
  125 + = Operators =
  126 + = =
  127 + ==========================================================================*/
  128 +
  129 +/***********************************************************************//**
  130 + * @brief Assignment operator
  131 + *
  132 + * @param[in] obs CTA on-off observation.
  133 + ***************************************************************************/
  134 +GCTAOnOffObservation& GCTAOnOffObservation::operator=(const GCTAOnOffObservation& obs)
  135 +{
  136 + // Execute only if object is not identical
  137 + if (this != &obs) {
  138 +
  139 + // Free members
  140 + free_members();
  141 +
  142 + // Initialise private members for clean destruction
  143 + init_members();
  144 +
  145 + // Copy members
  146 + copy_members(obs);
  147 +
  148 + } // endif: object was not identical
  149 +
  150 + // Return
  151 + return *this;
  152 +}
  153 +
  154 +
  155 +/*==========================================================================
  156 + = =
  157 + = Public methods =
  158 + = =
  159 + ==========================================================================*/
  160 +
  161 +/***********************************************************************//**
  162 + * @brief Clear instance
  163 + *
  164 + * This method properly resets the object to an initial state.
  165 + ***************************************************************************/
  166 +void GCTAOnOffObservation::clear(void)
  167 +{
  168 + // Free class members
  169 + free_members();
  170 +
  171 + // Initialise members
  172 + init_members();
  173 +
  174 + // Return
  175 + return;
  176 +}
  177 +
  178 +
  179 +/***********************************************************************//**
  180 + * @brief Clone instance
  181 + **************************************************************************/
  182 +GCTAOnOffObservation* GCTAOnOffObservation::clone(void) const
  183 +{
  184 + return new GCTAOnOffObservation(*this);
  185 +}
  186 +
  187 +
  188 +/***********************************************************************//**
  189 + * @brief read observation from an xml element
  190 + *
  191 + * @param[in] xml GXmlElement.
  192 + ***************************************************************************/
  193 +void GCTAOnOffObservation::read(const GXmlElement& xml)
  194 +{
  195 + // clean object
  196 + clear();
  197 +
  198 + // Extract instrument name
  199 + m_instrument = xml.attribute("instrument");
  200 +
  201 + // Determine number of parameter nodes in XML element
  202 + int npars = xml.elements("parameter");
  203 +
  204 + // Verify that XML element has exactly 6 parameters
  205 + if (xml.elements() != 6 || npars != 6) {
  206 + throw GException::xml_invalid_parnum(G_READ, xml,
  207 + "CTA observation requires exactly 6 parameters.");
  208 + }
  209 +
  210 + // Extract parameters
  211 + int npar[] = {0, 0, 0, 0, 0, 0};
  212 + for (int i = 0; i < npars; ++i) {
  213 +
  214 + // Get parameter element
  215 + const GXmlElement* par = xml.element("parameter", i);
  216 +
  217 + // Handle Pha_on
  218 + if (par->attribute("name") == "Pha_on") {
  219 +
  220 + // Read pha_on file name
  221 + std::string filename = par->attribute("file");
  222 +
  223 + // Load pha_on
  224 + m_on_spec.load(filename);
  225 +
  226 + // Increment parameter counter
  227 + npar[0]++;
  228 + }
  229 +
  230 + // Handle Pha_off
  231 + else if (par->attribute("name") == "Pha_off") {
  232 +
  233 + // Read countsmap file name
  234 + std::string filename = par->attribute("file");
  235 +
  236 + // Load pha_off
  237 + m_off_spec.load(filename);
  238 +
  239 + // Increment parameter counter
  240 + npar[1]++;
  241 + }
  242 +
  243 + // Handle on regions
  244 + else if(par->attribute("name") == "Regions_on") {
  245 +
  246 + // Get filename
  247 + std::string filename = par->attribute("file");
  248 +
  249 + // load on regions
  250 + m_on_regions.load(filename);
  251 +
  252 + // Increase number of parameters
  253 + npar[2]++;
  254 + }
  255 +
  256 + // Handle off regions
  257 + else if(par->attribute("name") == "Regions_off") {
  258 +
  259 + // Get filename
  260 + std::string filename = par->attribute("file");
  261 +
  262 + // load off regions
  263 + m_on_regions.load(filename);
  264 +
  265 + // Increase number of parameters
  266 + npar[3]++;
  267 + }
  268 +
  269 + // Handle Arf
  270 + else if(par->attribute("name") == "Arf") {
  271 +
  272 + // Get filename
  273 + std::string filename = par->attribute("file");
  274 +
  275 + // load arf
  276 + m_arf.load(filename);
  277 +
  278 + // Increase number of parameters
  279 + npar[4]++;
  280 + }
  281 +
  282 + // Handle Rmf
  283 + else if(par->attribute("name") == "Rmf") {
  284 +
  285 + // Get filename
  286 + std::string filename = par->attribute("file");
  287 +
  288 + // load Rmf
  289 + m_rmf.load(filename);
  290 +
  291 + // Increase number of parameters
  292 + npar[5]++;
  293 + }
  294 +
  295 + } // endfor: looped over all parameters
  296 +
  297 + // Verify that all parameters were found
  298 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1 || npar[4] != 1 || npar[5] != 1) {
  299 + throw GException::xml_invalid_parnames(G_READ, xml,
  300 + "Require \"Pha_on\" or \"Pha_off\" and \"Regions_on\""
  301 + ", \"Regions_off\",\"Rmf\" and \"Arf\" parameters.");
  302 + }
  303 +
  304 + // Return
  305 + return;
  306 +}
  307 +
  308 +
  309 +/***********************************************************************//**
  310 + * @brief write observation to an xml element
  311 + *
  312 + * @param[in] xml GXmlElement.
  313 + ***************************************************************************/
  314 +void GCTAOnOffObservation::write(GXmlElement& xml) const
  315 +{
  316 + // If XML element has 0 nodes then append 6 parameter nodes
  317 + if (xml.elements() == 0) {
  318 + xml.append(GXmlElement("parameter name=\"Pha_on\""));
  319 + xml.append(GXmlElement("parameter name=\"Pha_off\""));
  320 + xml.append(GXmlElement("parameter name=\"Regions_on\""));
  321 + xml.append(GXmlElement("parameter name=\"Regions_off\""));
  322 + xml.append(GXmlElement("parameter name=\"Arf\""));
  323 + xml.append(GXmlElement("parameter name=\"Rmf\""));
  324 + }
  325 +
  326 + // Verify that XML element has exactly 4 parameters
  327 + if (xml.elements() != 6 || xml.elements("parameter") != 6) {
  328 + throw GException::xml_invalid_parnum(G_WRITE, xml,
  329 + "CTAOnOffObservation requires exactly 6 parameters.");
  330 + }
  331 +
  332 + // Set or update parameter attributes
  333 + int npar[] = {0, 0, 0, 0, 0, 0};
  334 + for (int i = 0; i < 6; ++i) {
  335 +
  336 + // Get parameter element
  337 + GXmlElement* par = xml.element("parameter", i);
  338 +
  339 + // Handle on counts
  340 + if (par->attribute("name") == "Pha_on") {
  341 + par->attribute("file", m_on_spec.filename());
  342 + npar[0]++;
  343 + }
  344 +
  345 + // Handle off counts
  346 + else if (par->attribute("name") == "Pha_off") {
  347 + par->attribute("file", m_off_spec.filename());
  348 + npar[0]++;
  349 + }
  350 +
  351 + // Handle on regions
  352 + else if (par->attribute("name") == "Regions_on") {
  353 + par->attribute("file", m_on_regions.filename());
  354 + npar[0]++;
  355 + }
  356 +
  357 + // Handle off regions
  358 + else if (par->attribute("name") == "Regions_off") {
  359 + par->attribute("file", m_off_regions.filename());
  360 + npar[0]++;
  361 + }
  362 +
  363 + // Handle effective area
  364 + else if (par->attribute("name") == "Arf") {
  365 + par->attribute("file", m_arf.filename());
  366 + npar[0]++;
  367 + }
  368 +
  369 + // Handle energy resolution
  370 + else if (par->attribute("name") == "Rmf") {
  371 + par->attribute("file", m_rmf.filename());
  372 + npar[0]++;
  373 + }
  374 + }
  375 + // Verify that all required parameters are present
  376 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1 || npar[4] != 1 || npar[5] !=1) {
  377 + throw GException::xml_invalid_parnames(G_WRITE, xml,
  378 + "Require \"Pha_on\" or \"Pha_off\" and \"Regions_on\", \"Regions_off\""
  379 + " \"Arf\", \"Rmf\" parameters.");
  380 + }
  381 +
  382 + // Return
  383 + return;
  384 +}
  385 +
  386 +
  387 +/***********************************************************************//**
  388 + * @brief Fill events in the on and off spectra
  389 + *
  390 + * @param[in] obs CTA observation.
  391 + *
  392 + * @exception GException::invalid_value
  393 + * No CTA event list found in CTA observation.
  394 + ***************************************************************************/
  395 +void GCTAOnOffObservation::fill(const GCTAObservation& obs)
  396 +{
  397 + // Get CTA event list pointer
  398 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  399 + if (events == NULL) {
  400 + throw GException::invalid_value(G_FILL,
  401 + "GCTAObservation must provide an event list");
  402 + }
  403 +
  404 + // Loop over all events
  405 + for (int i = 0; i < events->size(); ++i) {
  406 +
  407 + // Get measured event direction
  408 + const GCTAEventAtom* atom = (*events)[i];
  409 + GSkyDir dir = atom->dir().dir();
  410 +
  411 + // Fill in spectrum according to region containment
  412 + if (m_on_regions.contains(dir)) {
  413 + m_on_spec.fill(atom->energy());
  414 + }
  415 + else if (m_off_regions.contains(dir)) {
  416 + m_off_spec.fill(atom->energy());
  417 + }
  418 +
  419 + } // endfor: looped over all events
  420 +
  421 + // Return
  422 + return;
  423 +}
  424 +
  425 +
  426 +/***********************************************************************//**
  427 + * @brief Compute response of an on/off observation
  428 + *
  429 + * @param[in] obs GCTAObservation.
  430 + * @param[in] etrue True energy boundaries.
  431 + ***************************************************************************/
  432 +void GCTAOnOffObservation::compute_response(const GCTAObservation& obs,
  433 + const GEbounds& etrue)
  434 +{
  435 + // Compute response components
  436 + compute_arf(obs);
  437 + compute_rmf(obs, etrue);
  438 +
  439 + // Return
  440 + return;
  441 +}
  442 +
  443 +/***********************************************************************//**
  444 + * @brief Compute arf of an on/off observation
  445 + *
  446 + * @param[in] obs GCTAObservation.
  447 + *
  448 + * @todo Implement GCTAResponse::npred usage.
  449 + ***************************************************************************/
  450 +void GCTAOnOffObservation::compute_arf(const GCTAObservation& obs)
  451 +{
  452 + // Set constant response parameters
  453 + const double theta = 0.0;
  454 + const double phi = 0.0;
  455 + const double zenith = 0.0;
  456 + const double azimuth = 0.0;
  457 +
  458 + // Get energy boundaries from on spectrum
  459 + GEbounds ereco = m_on_spec.ebounds();
  460 +
  461 + // Continue only if there are spectral bins
  462 + int nreco = ereco.size();
  463 + if (nreco > 0) {
  464 +
  465 + // Get CTA response pointer
  466 + GCTAResponse* response = obs.response();
  467 +
  468 + // Initialize ARF
  469 + m_arf = GArf(ereco);
  470 +
  471 + // Loop over reconstructed energies
  472 + for (int i = 0; i < nreco; ++i) {
  473 +
  474 + // Get mean energy of bin
  475 + double logenergy = ereco.elogmean(i).log10TeV();
  476 +
  477 + // Set specresp value
  478 + m_arf[i] = response->aeff(theta,
  479 + phi,
  480 + zenith,
  481 + azimuth,
  482 + logenergy);
  483 +
  484 + } // endfor: looped over reconstructed energies
  485 +
  486 + } // endif: there were energy bins
  487 +
  488 + // Return
  489 + return;
  490 +}
  491 +
  492 +
  493 +/***********************************************************************//**
  494 + * @brief compute rmf of an on/off observation
  495 + *
  496 + * @param[in] obs GCTAObservation.
  497 + * @param[in] etrue True energy boundaries.
  498 + ***************************************************************************/
  499 +void GCTAOnOffObservation::compute_rmf(const GCTAObservation& obs,
  500 + const GEbounds& etrue)
  501 +{
  502 + // Set constant response parameters
  503 + const double theta = 0.0;
  504 + const double phi = 0.0;
  505 + const double zenith = 0.0;
  506 + const double azimuth = 0.0;
  507 +
  508 + // Get reconstructed energy boundaries from on spectrum
  509 + GEbounds ereco = m_on_spec.ebounds();
  510 +
  511 + // Continue only if there are spectral bins
  512 + int ntrue = etrue.size();
  513 + int nreco = ereco.size();
  514 + if (ntrue > 0 && nreco > 0) {
  515 +
  516 + // Get CTA response pointer
  517 + GCTAResponse* response = obs.response();
  518 +
  519 + // Initialize RMF
  520 + m_rmf = GRmf(etrue, ereco);
  521 +
  522 + // Loop over reconstructed energy
  523 + for (int ireco = 0; ireco < nreco; ++ireco) {
  524 +
  525 + // Compute reconstructed energy
  526 + double eng_reco = ereco.elogmean(ireco).log10TeV();
  527 +
  528 + // Loop over true energy
  529 + for (int itrue = 0; itrue < ntrue; ++itrue) {
  530 +
  531 + // Compute true energy
  532 + double eng_true = etrue.elogmean(itrue).log10TeV();
  533 +
  534 + // Set RMF value
  535 + m_rmf(itrue, ireco) = response->edisp(eng_reco,
  536 + theta,
  537 + phi,
  538 + zenith,
  539 + azimuth,
  540 + eng_true);
  541 + } // endfor: looped over true energy
  542 + } // endfor: looped over reconstructed energy
  543 + } // endif: there were energy bins
  544 +
  545 + // Return
  546 + return;
  547 +}
  548 +
  549 +
  550 +/***********************************************************************//**
  551 + * @brief Print class information
  552 + ***************************************************************************/
  553 +std::string GCTAOnOffObservation::print(const GChatter& chatter) const
  554 +{
  555 + // Initialise result string
  556 + std::string result;
  557 +
  558 + // Continue only if chatter is not silent
  559 + if (chatter != SILENT) {
  560 +
  561 + // Append header
  562 + result.append("=== GCTAOnOffObservation ===");
  563 +
  564 + // Append parameters
  565 + result.append("\n"+gammalib::parformat("Name")+m_name);
  566 + result.append("\n"+gammalib::parformat("Identifier")+m_id);
  567 +
  568 + // Append spectra, ARF and RMF
  569 + result.append("\n"+m_on_spec.print());
  570 + result.append("\n"+m_off_spec.print());
  571 + result.append("\n"+m_arf.print());
  572 + result.append("\n"+m_rmf.print());
  573 +
  574 + // Append regions
  575 + result.append("\n"+m_on_regions.print());
  576 + result.append("\n"+m_off_regions.print());
  577 + }
  578 +
  579 + // Return result
  580 + return result;
  581 +}
  582 +
  583 +
  584 +/*==========================================================================
  585 + = =
  586 + = Private methods =
  587 + = =
  588 + ==========================================================================*/
  589 +
  590 +/***********************************************************************//**
  591 + * @brief Initialise class members
  592 + ***************************************************************************/
  593 +void GCTAOnOffObservation::init_members(void)
  594 +{
  595 + // Initialise members
  596 + m_name.clear();
  597 + m_id.clear();
  598 + m_instrument = "CTA";
  599 + m_on_spec.clear();
  600 + m_off_spec.clear();
  601 + m_arf.clear();
  602 + m_rmf.clear();
  603 + m_on_regions.clear();
  604 + m_off_regions.clear();
  605 +
  606 + // Return
  607 + return;
  608 +}
  609 +
  610 +
  611 +/***********************************************************************//**
  612 + * @brief Copy class members
  613 + *
  614 + * @param[in] obs CTA on-off observation.
  615 + ***************************************************************************/
  616 +void GCTAOnOffObservation::copy_members(const GCTAOnOffObservation& obs)
  617 +{
  618 + // Copy attributes
  619 + m_name = obs.m_name;
  620 + m_id = obs.m_id;
  621 + m_instrument = obs.m_instrument;
  622 + m_on_spec = obs.m_on_spec;
  623 + m_off_spec = obs.m_off_spec;
  624 + m_arf = obs.m_arf;
  625 + m_rmf = obs.m_rmf;
  626 + m_on_regions = obs.m_on_regions;
  627 + m_off_regions = obs.m_off_regions;
  628 +
  629 + // Return
  630 + return;
  631 +}
  632 +
  633 +
  634 +
  635 +/***********************************************************************//**
  636 + * @brief Delete class members
  637 + ***************************************************************************/
  638 +void GCTAOnOffObservation::free_members(void)
  639 +{
  640 + // Return
  641 + return;
  642 +}
... ...
inst/cta/src/GCTAOnOffObservations.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GCTAOnOffObservations.cpp - ON/OFF Observation container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservations.cpp
  23 + * @brief ON/OFF Observations container class implementation
  24 + * @author Pierrick Martin
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GException.hpp"
  33 +#include "GCTAOnOffObservations.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_AT "GCTAOnOffObservations::at(int&)"
  37 +#define G_SET "GCTAOnOffObservations::set(int&, GCTAOnOffObservation&)"
  38 +#define G_APPEND "GCTAOnOffObservations::append(GCTAOnOffObservation&)"
  39 +#define G_INSERT "GCTAOnOffObservations::insert(int&, GCTAOnOffObservation&)"
  40 +#define G_REMOVE "GCTAOnOffObservations::remove(int&)"
  41 +#define G_EXTEND "GCTAOnOffObservations::extend(GCTAOnOffObservations&)"
  42 +#define G_READ "GCTAOnOffObservations::read(GXml&)"
  43 +
  44 +/* __ Macros _____________________________________________________________ */
  45 +
  46 +/* __ Coding definitions _________________________________________________ */
  47 +
  48 +/* __ Debug definitions __________________________________________________ */
  49 +
  50 +
  51 +/*==========================================================================
  52 + = =
  53 + = Constructors/destructors =
  54 + = =
  55 + ==========================================================================*/
  56 +
  57 +/***********************************************************************//**
  58 + * @brief Void constructor
  59 + ***************************************************************************/
  60 +GCTAOnOffObservations::GCTAOnOffObservations(void)
  61 +{
  62 + // Initialise members
  63 + init_members();
  64 +
  65 + // Return
  66 + return;
  67 +}
  68 +
  69 +
  70 +/***********************************************************************//**
  71 + * @brief Copy constructor
  72 + *
  73 + * @param obs Observation container.
  74 + ***************************************************************************/
  75 +GCTAOnOffObservations::GCTAOnOffObservations(const GCTAOnOffObservations& obs)
  76 +{
  77 + // Initialise members
  78 + init_members();
  79 +
  80 + // Copy members
  81 + copy_members(obs);
  82 +
  83 + // Return
  84 + return;
  85 +}
  86 +
  87 +
  88 +/***********************************************************************//**
  89 + * @brief Load constructor
  90 + *
  91 + * @param[in] filename XML filename.
  92 + *
  93 + * Construct the observation container by loading all relevant information
  94 + * from a XML file. Please refer to the read() method for more information
  95 + * about the structure of the XML file.
  96 + ***************************************************************************/
  97 +GCTAOnOffObservations::GCTAOnOffObservations(const std::string& filename)
  98 +{
  99 + // Initialise members
  100 + init_members();
  101 +
  102 + // Load XML file
  103 + load(filename);
  104 +
  105 + // Return
  106 + return;
  107 +}
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Destructor
  112 + ***************************************************************************/
  113 +GCTAOnOffObservations::~GCTAOnOffObservations(void)
  114 +{
  115 + // Free members
  116 + free_members();
  117 +
  118 + // Return
  119 + return;
  120 +}
  121 +
  122 +
  123 +/*==========================================================================
  124 + = =
  125 + = Operators =
  126 + = =
  127 + ==========================================================================*/
  128 +
  129 +/***********************************************************************//**
  130 + * @brief Assignment operator
  131 + *
  132 + * @param[in] obs Observation container.
  133 + * @return Observation container.
  134 + ***************************************************************************/
  135 +GCTAOnOffObservations& GCTAOnOffObservations::operator=(const GCTAOnOffObservations& obs)
  136 +{
  137 + // Execute only if object is not identical
  138 + if (this != &obs) {
  139 +
  140 + // Free members
  141 + free_members();
  142 +
  143 + // Initialise members
  144 + init_members();
  145 +
  146 + // Copy members
  147 + copy_members(obs);
  148 +
  149 + } // endif: object was not identical
  150 +
  151 + // Return this object
  152 + return *this;
  153 +}
  154 +
  155 +
  156 +/*==========================================================================
  157 + = =
  158 + = Public methods =
  159 + = =
  160 + ==========================================================================*/
  161 +
  162 +/***********************************************************************//**
  163 + * @brief Clear observations
  164 + ***************************************************************************/
  165 +void GCTAOnOffObservations::clear(void)
  166 +{
  167 + // Free members
  168 + free_members();
  169 +
  170 + // Initialise members
  171 + init_members();
  172 +
  173 + // Return
  174 + return;
  175 +}
  176 +
  177 +
  178 +/***********************************************************************//**
  179 + * @brief Clone observations
  180 + *
  181 + * @return Pointer to deep copy of observation container.
  182 + ***************************************************************************/
  183 +GCTAOnOffObservations* GCTAOnOffObservations::clone(void) const
  184 +{
  185 + // Clone observations
  186 + return new GCTAOnOffObservations(*this);
  187 +}
  188 +
  189 +
  190 +/***********************************************************************//**
  191 + * @brief Return pointer to observation
  192 + *
  193 + * @param[in] index Observation index [0,...,size()-1].
  194 + * @return Observation.
  195 + *
  196 + * @exception GException::out_of_range
  197 + * Operation index is out of range.
  198 + *
  199 + * Returns a pointer to the observation with specified @p index.
  200 + ***************************************************************************/
  201 +GCTAOnOffObservation* GCTAOnOffObservations::at(const int& index)
  202 +{
  203 + // Raise exception if index is out of range
  204 + if (index < 0 || index >= size()) {
  205 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  206 + }
  207 +
  208 + // Return pointer
  209 + return m_obs[index];
  210 +}
  211 +
  212 +
  213 +/***********************************************************************//**
  214 + * @brief Return pointer to observation (const version)
  215 + *
  216 + * @param[in] index Observation index [0,...,size()-1].
  217 + *
  218 + * @exception GException::out_of_range
  219 + * Operation index is out of range.
  220 + *
  221 + * Returns a const pointer to the observation with specified @p index.
  222 + ***************************************************************************/
  223 +const GCTAOnOffObservation* GCTAOnOffObservations::at(const int& index) const
  224 +{
  225 + // Raise exception if index is out of range
  226 + if (index < 0 || index >= size()) {
  227 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  228 + }
  229 +
  230 + // Return pointer
  231 + return m_obs[index];
  232 +}
  233 +
  234 +
  235 +/***********************************************************************//**
  236 + * @brief Set observation in container
  237 + *
  238 + * @param[in] index Observation index [0,...,size()-1].
  239 + * @param[in] obs Observation.
  240 + * @return Pointer to deep copy of observation.
  241 + *
  242 + * @exception GException::out_of_range
  243 + * Observation index is out of range.
  244 + * @exception GException::invalid_value
  245 + * Observation with same instrument and identifier already
  246 + * exists in container.
  247 + *
  248 + * Set a deep copy and observation @p obs at the specified @p index in the
  249 + * container.
  250 + ***************************************************************************/
  251 +GCTAOnOffObservation* GCTAOnOffObservations::set(const int& index, const GCTAOnOffObservation& obs)
  252 +{
  253 + // Compile option: raise exception if index is out of range
  254 + #if defined(G_RANGE_CHECK)
  255 + if (index < 0 || index >= size()) {
  256 + throw GException::out_of_range(G_SET, index, 0, size()-1);
  257 + }
  258 + #endif
  259 +
  260 + // Raise an exception if an observation with specified instrument and
  261 + // identifier already exists
  262 + int inx = get_index(obs.instrument(), obs.id());
  263 + if (inx != -1 && inx != index) {
  264 + std::string msg =
  265 + "Attempt to set \""+obs.instrument()+"\" observation with"
  266 + " identifier \""+obs.id()+"\" in observation container at"
  267 + " index "+gammalib::str(index)+", but an observation with the"
  268 + " same attributes exists already at index "+gammalib::str(inx)+
  269 + " in the container.\n"
  270 + "Every observation for a given instrument in the observation"
  271 + " container needs a unique identifier.";
  272 + throw GException::invalid_value(G_SET, msg);
  273 + }
  274 +
  275 + // Delete existing observation
  276 + if (m_obs[index] != NULL) delete m_obs[index];
  277 +
  278 + // Store pointer to a deep copy of the observation
  279 + m_obs[index] = obs.clone();
  280 +
  281 + // Return pointer
  282 + return m_obs[index];
  283 +}
  284 +
  285 +
  286 +/***********************************************************************//**
  287 + * @brief Append observation to container
  288 + *
  289 + * @param[in] obs Observation.
  290 + * @return Pointer to deep copy of observation.
  291 + *
  292 + * @exception GException::invalid_value
  293 + * Observation with same instrument and identifier already
  294 + * exists in container.
  295 + *
  296 + * Appends a deep copy of an observation to the container.
  297 + ***************************************************************************/
  298 +GCTAOnOffObservation* GCTAOnOffObservations::append(const GCTAOnOffObservation& obs)
  299 +{
  300 + // Raise an exception if an observation with specified instrument and
  301 + // identifier already exists
  302 + int inx = get_index(obs.instrument(), obs.id());
  303 + if (inx != -1) {
  304 + std::string msg =
  305 + "Attempt to append \""+obs.instrument()+"\" observation with"
  306 + " identifier \""+obs.id()+"\" to observation container, but an"
  307 + " observation with the same attributes exists already at"
  308 + " index "+gammalib::str(inx)+" in the container.\n"
  309 + "Every observation for a given instrument in the observation"
  310 + " container needs a unique identifier.";
  311 + throw GException::invalid_value(G_APPEND, msg);
  312 + }
  313 +
  314 + // Clone observation
  315 + GCTAOnOffObservation* ptr = obs.clone();
  316 +
  317 + // Append to list
  318 + m_obs.push_back(ptr);
  319 +
  320 + // Return pointer
  321 + return ptr;
  322 +}
  323 +
  324 +
  325 +/***********************************************************************//**
  326 + * @brief Insert observation into container
  327 + *
  328 + * @param[in] index Observation index [0,...,size()-1].
  329 + * @param[in] obs Observation.
  330 + * @return Pointer to deep copy of observation.
  331 + *
  332 + * @exception GException::out_of_range
  333 + * Observation index is out of range.
  334 + * @exception GException::invalid_value
  335 + * Observation with same instrument and identifier already
  336 + * exists in container.
  337 + *
  338 + * Inserts a deep copy of an observation into the container before the
  339 + * observation with the specified @p index.
  340 + ***************************************************************************/
  341 +GCTAOnOffObservation* GCTAOnOffObservations::insert(const int& index, const GCTAOnOffObservation& obs)
  342 +{
  343 + // Compile option: raise exception if index is out of range
  344 + #if defined(G_RANGE_CHECK)
  345 + if (isempty()) {
  346 + if (index > 0) {
  347 + throw GException::out_of_range(G_INSERT, index, 0, size()-1);
  348 + }
  349 + }
  350 + else {
  351 + if (index < 0 || index >= size()) {
  352 + throw GException::out_of_range(G_INSERT, index, 0, size()-1);
  353 + }
  354 + }
  355 + #endif
  356 +
  357 + // Raise an exception if an observation with specified instrument and
  358 + // identifier already exists
  359 + int inx = get_index(obs.instrument(), obs.id());
  360 + if (inx != -1 && inx != index) {
  361 + std::string msg =
  362 + "Attempt to insert \""+obs.instrument()+"\" observation with"
  363 + " identifier \""+obs.id()+"\" in observation container before"
  364 + " index "+gammalib::str(index)+", but an observation with the"
  365 + " same attributes exists already at index "+gammalib::str(inx)+
  366 + " in the container.\n"
  367 + "Every observation for a given instrument in the observation"
  368 + " container needs a unique identifier.";
  369 + throw GException::invalid_value(G_INSERT, msg);
  370 + }
  371 +
  372 + // Clone observation
  373 + GCTAOnOffObservation* ptr = obs.clone();
  374 +
  375 + // Clone observation and insert into list
  376 + m_obs.insert(m_obs.begin()+index, ptr);
  377 +
  378 + // Return pointer
  379 + return ptr;
  380 +}
  381 +
  382 +
  383 +/***********************************************************************//**
  384 + * @brief Remove observation from container
  385 + *
  386 + * @param[in] index Observation index [0,...,size()-1].
  387 + *
  388 + * @exception GException::out_of_range
  389 + * Observation index is out of range.
  390 + *
  391 + * Removes observation of specified @p index from the container.
  392 + ***************************************************************************/
  393 +void GCTAOnOffObservations::remove(const int& index)
  394 +{
  395 + // Compile option: If index is outside boundary then raise exception
  396 + #if defined(G_RANGE_CHECK)
  397 + if (index < 0 || index >= size()) {
  398 + throw GException::out_of_range(G_REMOVE, index, 0, size()-1);
  399 + }
  400 + #endif
  401 +
  402 + // Erase observation component from container
  403 + m_obs.erase(m_obs.begin() + index);
  404 +
  405 + // Return
  406 + return;
  407 +}
  408 +
  409 +
  410 +/***********************************************************************//**
  411 + * @brief Append observations from observation container
  412 + *
  413 + * @param[in] obs Observations.
  414 + *
  415 + * @exception GException::invalid_value
  416 + * Observation with same instrument and identifier already
  417 + * exists in container.
  418 + *
  419 + * Appends deep copies of observations to the container.
  420 + ***************************************************************************/
  421 +void GCTAOnOffObservations::extend(const GCTAOnOffObservations& obs)
  422 +{
  423 + // Do nothing if observation container is empty
  424 + if (!obs.isempty()) {
  425 +
  426 + // Get size. Note that we extract the size first to avoid an
  427 + // endless loop that arises when a container is appended to
  428 + // itself.
  429 + int num = obs.size();
  430 +
  431 + // Reserve enough space
  432 + reserve(size() + num);
  433 +
  434 + // Loop over all observations, clone them and append them to the
  435 + // list
  436 + for (int i = 0; i < num; ++i) {
  437 +
  438 + // Raise an exception if an observation with specified
  439 + // instrument and identifier already exists
  440 + int inx = get_index(obs[i]->instrument(), obs[i]->id());
  441 + if (inx != -1) {
  442 + std::string msg =
  443 + "Attempt to append \""+obs[i]->instrument()+"\""
  444 + " observation with identifier \""+obs[i]->id()+"\""
  445 + " to observation container, but an observation with the"
  446 + " same attributes exists already at index "+
  447 + gammalib::str(inx)+" in the container.\n"
  448 + "Every observation for a given instrument in the"
  449 + " observation container needs a unique identifier.";
  450 + throw GException::invalid_value(G_EXTEND, msg);
  451 + }
  452 +
  453 + // Append observation to container
  454 + m_obs.push_back(obs[i]->clone());
  455 + }
  456 +
  457 + } // endif: observation container was not empty
  458 +
  459 + // Return
  460 + return;
  461 +}
  462 +
  463 +
  464 +/***********************************************************************//**
  465 + * @brief Signals if observation exists
  466 + *
  467 + * @param[in] instrument Instrument.
  468 + * @param[in] id Observation identifier.
  469 + * @return True if observation with specified @p instrument and identifier
  470 + * @p id exists.
  471 + *
  472 + * Searches all observations for a match with the specified @p instrument
  473 + * and @p identifier. If the specified attributes have been found, true is
  474 + * returned.
  475 + ***************************************************************************/
  476 +bool GCTAOnOffObservations::contains(const std::string& instrument,
  477 + const std::string& id) const
  478 +{
  479 + // Get observation index
  480 + int index = get_index(instrument, id);
  481 +
  482 + // Return
  483 + return (index != -1);
  484 +}
  485 +
  486 +
  487 +/***********************************************************************//**
  488 + * @brief Load observations from XML file
  489 + *
  490 + * @param[in] filename XML filename.
  491 + *
  492 + * Loads observation from a XML file into the container. Please refer to the
  493 + * read() method for more information about the structure of the XML file.
  494 + ***************************************************************************/
  495 +void GCTAOnOffObservations::load(const std::string& filename)
  496 +{
  497 + // Clear any existing observations
  498 + clear();
  499 +
  500 + // Load XML document
  501 + GXml xml(filename);
  502 +
  503 + // Read observations from XML document
  504 + read(xml);
  505 +
  506 + // Return
  507 + return;
  508 +}
  509 +
  510 +
  511 +/***********************************************************************//**
  512 + * @brief Save observations into XML file
  513 + *
  514 + * @param[in] filename XML filename.
  515 + *
  516 + * Saves observations into a XML file. Please refer to the read() method for
  517 + * more information about the structure of the XML file.
  518 + ***************************************************************************/
  519 +void GCTAOnOffObservations::save(const std::string& filename) const
  520 +{
  521 + // Declare empty XML document
  522 + GXml xml;
  523 +
  524 + // Write observations into XML file
  525 + write(xml);
  526 +
  527 + // Save XML document
  528 + xml.save(filename);
  529 +
  530 + // Return
  531 + return;
  532 +}
  533 +
  534 +
  535 +/***********************************************************************//**
  536 + * @brief Read observations from XML document
  537 + *
  538 + * @param[in] xml XML document.
  539 + *
  540 + * @exception GException::invalid_instrument
  541 + * Invalid instrument encountered in XML file.
  542 + *
  543 + * Reads observations from the first observation list that is found in the
  544 + * XML document. The decoding of the instrument specific observation
  545 + * definition is done within the observation's GCTAOnOffObservation::read() method.
  546 + * The following file structure is expected:
  547 + *
  548 + * <observation_list title="observation library">
  549 + * <observation name="..." id="..." instrument="...">
  550 + * ...
  551 + * </observation>
  552 + * <observation name="..." id="..." instrument="...">
  553 + * ...
  554 + * </observation>
  555 + * ...
  556 + * </observation_list>
  557 + *
  558 + * The @p name and @p id attributes allow for a unique identification of an
  559 + * observation within the observation container. The @p instrument
  560 + * attributes specifies the instrument to which the observation applies.
  561 + * At the moment there is only CTA but more instruments may be implemented
  562 + * later and require the need for a registry (see GObservations.cpp).
  563 + *
  564 + * The structure within the @p observation tag is defined by the instrument
  565 + * specific GCTAOnOffObservation class.
  566 + *
  567 + ***************************************************************************/
  568 +void GCTAOnOffObservations::read(const GXml& xml)
  569 +{
  570 + // Get pointer on observation library
  571 + const GXmlElement* lib = xml.element("observation_list", 0);
  572 +
  573 + // Loop over all observations
  574 + int n = lib->elements("observation");
  575 + for (int i = 0; i < n; ++i) {
  576 +
  577 + // Get pointer on observation
  578 + const GXmlElement* obs = lib->element("observation", i);
  579 +
  580 + // Get attributes
  581 + std::string name = obs->attribute("name");
  582 + std::string id = obs->attribute("id");
  583 + std::string instrument = obs->attribute("instrument");
  584 +
  585 + // Allocate observation
  586 + // (only CTA at the moment, implement registry if more)
  587 + GCTAOnOffObservation* ptr = new GCTAOnOffObservation;
  588 +
  589 + // If observation is valid then read its definition from XML file
  590 + if (ptr != NULL) {
  591 +
  592 + // Read definition
  593 + ptr->read(*obs);
  594 +
  595 + // Set attributes
  596 + ptr->name(name);
  597 + ptr->id(id);
  598 +
  599 + } // endif: observation was valid
  600 +
  601 + // ... otherwise throw an exception
  602 + else {
  603 + std::string msg =
  604 + "Problem allocating GCTAOnOffObservation object.";
  605 + throw GException::invalid_instrument(G_READ, msg);
  606 + }
  607 +
  608 + // Append observation to container
  609 + append(*ptr);
  610 +
  611 + // Free observation (the append method made a deep copy)
  612 + delete ptr;
  613 +
  614 + } // endfor: looped over all observations
  615 +
  616 + // Return
  617 + return;
  618 +}
  619 +
  620 +
  621 +/***********************************************************************//**
  622 + * @brief Write observations into XML document
  623 + *
  624 + * @param[in] xml XML document.
  625 + *
  626 + * Write observations into the first observation library that is found in the
  627 + * XML document. In case that no observation library exists, one is added to
  628 + * the document. Please refer to the read() method for more information
  629 + * about the structure of the XML document.
  630 + ***************************************************************************/
  631 +void GCTAOnOffObservations::write(GXml& xml) const
  632 +{
  633 + // If there is no observation library then append one
  634 + if (xml.elements("observation_list") == 0) {
  635 + xml.append(GXmlElement("observation_list title=\"observation list\""));
  636 + }
  637 +
  638 + // Get pointer on observation library
  639 + GXmlElement* lib = xml.element("observation_list", 0);
  640 +
  641 + // Write all observations into library
  642 + for (int i = 0; i < size(); ++i) {
  643 +
  644 + // Initialise pointer on observation
  645 + GXmlElement* obs = NULL;
  646 +
  647 + // Search corresponding observation
  648 + int n = xml.elements("observation");
  649 + for (int k = 0; k < n; ++k) {
  650 + GXmlElement* element = xml.element("observation", k);
  651 + if (element->attribute("name") == m_obs[i]->name() &&
  652 + element->attribute("id") == m_obs[i]->id() &&
  653 + element->attribute("instrument") == m_obs[i]->instrument()) {
  654 + obs = element;
  655 + break;
  656 + }
  657 + }
  658 +
  659 + // If no observation with corresponding name, ID and instrument was
  660 + // found then append one now
  661 + if (obs == NULL) {
  662 + obs = lib->append("observation");
  663 + obs->attribute("name", m_obs[i]->name());
  664 + obs->attribute("id", m_obs[i]->id());
  665 + obs->attribute("instrument", m_obs[i]->instrument());
  666 + }
  667 +
  668 + // Write now observation
  669 + m_obs[i]->write(*obs);
  670 +
  671 + } // endfor: looped over all observaitons
  672 +
  673 + // Return
  674 + return;
  675 +}
  676 +
  677 +
  678 +/***********************************************************************//**
  679 + * @brief Load models from XML file
  680 + *
  681 + * @param[in] filename XML filename.
  682 + *
  683 + * Loads the models from a XML file. Please refer to the GModels::read()
  684 + * method for more information about the expected structure of the XML
  685 + * file.
  686 + ***************************************************************************/
  687 +void GCTAOnOffObservations::models(const std::string& filename)
  688 +{
  689 + // Load models
  690 + m_models.load(filename);
  691 +
  692 + // Return
  693 + return;
  694 +}
  695 +
  696 +
  697 +/***********************************************************************//**
  698 + * @brief Print observation list information
  699 + *
  700 + * @param[in] chatter Chattiness (defaults to NORMAL).
  701 + * @return String containing observation list information
  702 + ***************************************************************************/
  703 +std::string GCTAOnOffObservations::print(const GChatter& chatter) const
  704 +{
  705 + // Initialise result string
  706 + std::string result;
  707 +
  708 + // Continue only if chatter is not silent
  709 + if (chatter != SILENT) {
  710 +
  711 + // Append header
  712 + result.append("=== GCTAOnOffObservations ===");
  713 +
  714 + // Append information
  715 + result.append("\n"+gammalib::parformat("Number of observations"));
  716 + result.append(gammalib::str(size()));
  717 +
  718 + // NORMAL: Append observations
  719 + if (chatter >= NORMAL) {
  720 + for (int i = 0; i < size(); ++i) {
  721 + result.append("\n");
  722 + result.append((*this)[i]->print());
  723 + }
  724 + } // endif: chatter was normal
  725 +
  726 + // EXPLICIT: Append models
  727 + if (chatter >= EXPLICIT) {
  728 + result.append("\n"+m_models.print());
  729 + } // endif: chatter was explicit
  730 +
  731 + } // endif: chatter was not silent
  732 +
  733 + // Return result
  734 + return result;
  735 +}
  736 +
  737 +
  738 +/*==========================================================================
  739 + = =
  740 + = Private methods =
  741 + = =
  742 + ==========================================================================*/
  743 +
  744 +/***********************************************************************//**
  745 + * @brief Initialise class members
  746 + ***************************************************************************/
  747 +void GCTAOnOffObservations::init_members(void)
  748 +{
  749 + // Initialise members
  750 + m_obs.clear();
  751 + m_models.clear();
  752 +
  753 + // Return
  754 + return;
  755 +}
  756 +
  757 +
  758 +/***********************************************************************//**
  759 + * @brief Copy class members
  760 + *
  761 + * @param[in] obs Observation container.
  762 + *
  763 + * Copies all class members. Deep copies are created for the observations,
  764 + * which is what is needed to have a real physical copy of the members.
  765 + ***************************************************************************/
  766 +void GCTAOnOffObservations::copy_members(const GCTAOnOffObservations& obs)
  767 +{
  768 + // Copy attributes
  769 + m_models = obs.m_models;
  770 +
  771 + // Copy observations
  772 + m_obs.clear();
  773 + for (int i = 0; i < obs.m_obs.size(); ++i) {
  774 + m_obs.push_back((obs.m_obs[i]->clone()));
  775 + }
  776 +
  777 + // Return
  778 + return;
  779 +}
  780 +
  781 +
  782 +/***********************************************************************//**
  783 + * @brief Delete class members
  784 + *
  785 + * Deallocates all observations. Since container classes that hold pointers
  786 + * need to handle the proper deallocation of memory, we loop here over all
  787 + * pointers and make sure that we deallocate all observations in the
  788 + * container.
  789 + ***************************************************************************/
  790 +void GCTAOnOffObservations::free_members(void)
  791 +{
  792 + // Free observations
  793 + for (int i = 0; i < m_obs.size(); ++i) {
  794 + delete m_obs[i];
  795 + m_obs[i] = NULL;
  796 + }
  797 +
  798 + // Return
  799 + return;
  800 +}
  801 +
  802 +
  803 +/***********************************************************************//**
  804 + * @brief Return observation index by instrument and identifier
  805 + *
  806 + * @param[in] instrument Instrument.
  807 + * @param[in] id Observation identifier.
  808 + * @return Observation index (-1 of instrument and observation identifier
  809 + * has not been found)
  810 + *
  811 + * Returns observation index based on the specified @p instrument and
  812 + * observation identifier @p id. If no observation with the specified
  813 + * attributes has been found, the method returns -1.
  814 + ***************************************************************************/
  815 +int GCTAOnOffObservations::get_index(const std::string& instrument,
  816 + const std::string& id) const
  817 +{
  818 + // Initialise index
  819 + int index = -1;
  820 +
  821 + // Search observation with specified instrument and id
  822 + for (int i = 0; i < size(); ++i) {
  823 + if ((m_obs[i]->instrument() == instrument) &&
  824 + (m_obs[i]->id() == id)) {
  825 + index = i;
  826 + break;
  827 + }
  828 + }
  829 +
  830 + // Return index
  831 + return index;
  832 +}
... ...
pyext/GArf.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GArf.i - XSPEC Auxiliary Response File class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GArf.i
  23 + * @brief XSPEC Auxiliary Response File class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GArf.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GArf
  35 + *
  36 + * @brief Auxiliary Response File class
  37 + ***************************************************************************/
  38 +class GArf : public GBase {
  39 +public:
  40 + // Constructors and destructors
  41 + GArf(void);
  42 + explicit GArf(const std::string& filename);
  43 + explicit GArf(const GEbounds& ebds);
  44 + GArf(const GArf& arf);
  45 + virtual ~GArf(void);
  46 +
  47 + // Methods
  48 + void clear(void);
  49 + GArf* clone(void) const;
  50 + int size(void) const;
  51 + double& at(const int& index);
  52 + const GEbounds& ebounds(void) const;
  53 + void load(const std::string& filename);
  54 + void save(const std::string& filename, const bool& clobber = false) const;
  55 + void read(const GFitsTable* hdu);
  56 + void write(GFits& fits) const;
  57 +};
  58 +
  59 +
  60 +/***********************************************************************//**
  61 + * @brief GArf class extension
  62 + ***************************************************************************/
  63 +%extend GArf {
  64 + double __getitem__(const int& index) {
  65 + if (index >= 0 && index < self->size()) {
  66 + return (*self)[index];
  67 + }
  68 + else {
  69 + throw GException::out_of_range("__getitem__(int)", index, self->size());
  70 + }
  71 + }
  72 + void __setitem__(const int& index, const double& value) {
  73 + if (index>=0 && index < self->size()) {
  74 + (*self)[index] = value;
  75 + return;
  76 + }
  77 + else {
  78 + throw GException::out_of_range("__setitem__(int)", index, self->size());
  79 + }
  80 + }
  81 + int __len__() {
  82 + return (self->size());
  83 + }
  84 + GArf copy() {
  85 + return (*self);
  86 + }
  87 +};
... ...
pyext/GPha.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GPha.i - XSPEC Pulse Height Analyzer class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GPha.i
  23 + * @brief XSPEC Pulse Height Analyzer class
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GPha.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GPha
  35 + *
  36 + * @brief Pulse Height Analyzer class
  37 + ***************************************************************************/
  38 +class GPha : public GBase {
  39 +public:
  40 + // Constructors and destructors
  41 + GPha(void);
  42 + explicit GPha(const std::string& filename);
  43 + explicit GPha(const GEbounds& ebds);
  44 + GPha(const GPha& pha);
  45 + virtual ~GPha(void);
  46 +
  47 + // Methods
  48 + void clear(void);
  49 + GPha* clone(void) const;
  50 + int size(void) const;
  51 + double& at(const int& index);
  52 + const GEbounds& ebounds(void) const;
  53 + double counts(void) const;
  54 + const double& underflow(void) const;
  55 + const double& overflow(void) const;
  56 + const double& outflow(void) const;
  57 + void fill(const GEnergy& energy, const double& value = 1.0);
  58 + void load(const std::string& filename);
  59 + void save(const std::string& filename, const bool& clobber = false) const;
  60 + void read(const GFitsTable* hdu);
  61 + void write(GFits& fits) const;
  62 +};
  63 +
  64 +
  65 +/***********************************************************************//**
  66 + * @brief GPha class extension
  67 + ***************************************************************************/
  68 +%extend GPha {
  69 + double __getitem__(const int& index) {
  70 + if (index >= 0 && index < self->size()) {
  71 + return (*self)[index];
  72 + }
  73 + else {
  74 + throw GException::out_of_range("__getitem__(int)", index, self->size());
  75 + }
  76 + }
  77 + void __setitem__(const int& index, const double& value) {
  78 + if (index>=0 && index < self->size()) {
  79 + (*self)[index] = value;
  80 + return;
  81 + }
  82 + else {
  83 + throw GException::out_of_range("__setitem__(int)", index, self->size());
  84 + }
  85 + }
  86 + int __len__() {
  87 + return (self->size());
  88 + }
  89 + GPha copy() {
  90 + return (*self);
  91 + }
  92 +};
... ...
pyext/GRmf.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GRmf.i - XSPEC Redistribution Matrix File class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GRmf.i
  23 + * @brief XSPEC Redistribution Matrix File class
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GRmf.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +/* __ Includes ___________________________________________________________ */
  33 +%include "GTypemaps.i"
  34 +
  35 +
  36 +/***********************************************************************//**
  37 + * @class GRmf
  38 + *
  39 + * @brief Redistribution Matrix File class
  40 + ***************************************************************************/
  41 +class GRmf : public GBase {
  42 +public:
  43 + // Constructors and destructors
  44 + GRmf(void);
  45 + explicit GRmf(const std::string& filename);
  46 + explicit GRmf(const GEbounds& etrue, const GEbounds& emeasured);
  47 + GRmf(const GRmf& rmf);
  48 + virtual ~GRmf(void);
  49 +
  50 + // Methods
  51 + void clear(void);
  52 + GRmf* clone(void) const;
  53 + int size(void) const;
  54 + int ntrue(void) const;
  55 + int nmeasured(void) const;
  56 + double& at(const int& itrue, const int& imeasured);
  57 + const GEbounds& etrue(void) const;
  58 + const GEbounds& emeasured(void) const;
  59 + void load(const std::string& filename);
  60 + void save(const std::string& filename, const bool& clobber = false) const;
  61 + void read(const GFitsTable* hdu);
  62 + void write(GFits& fits) const;
  63 +};
  64 +
  65 +
  66 +/***********************************************************************//**
  67 + * @brief GRmf class extension
  68 + ***************************************************************************/
  69 +%extend GRmf {
  70 + double __getitem__(int GTuple[2]) {
  71 + return self->at(GTuple[0], GTuple[1]);
  72 + }
  73 + void __setitem__(int GTuple[2], const double& value) {
  74 + self->at(GTuple[0], GTuple[1]) = value;
  75 + }
  76 + int __len__() {
  77 + return (self->size());
  78 + }
  79 + GRmf copy() {
  80 + return (*self);
  81 + }
  82 +};
... ...
pyext/GSkyRegion.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegion.i - Sky region class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GSkyRegion.i
  23 + * @brief Sky region class SWIG file.
  24 + * @author Michael Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GSkyRegion.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GSkyRegion
  35 + *
  36 + * @brief GSkyRegion base abstract class
  37 + ***************************************************************************/
  38 +class GSkyRegion : public GBase {
  39 +public:
  40 + // Constructors and destructors
  41 + GSkyRegion(void);
  42 + GSkyRegion(const GSkyRegion& region);
  43 + virtual ~GSkyRegion(void);
  44 +
  45 + // Pure virtual methods
  46 + virtual void clear(void) = 0;
  47 + virtual GSkyRegion* clone(void) const = 0;
  48 + virtual void read(const std::string regstring) = 0;
  49 + virtual std::string write(void) const = 0;
  50 + virtual bool contains(const GSkyDir& dir) const = 0;
  51 + virtual bool overlaps(const GSkyRegion& reg) const = 0;
  52 + virtual bool contains(const GSkyRegion& reg) const = 0;
  53 +
  54 + // Implemented methods
  55 + std::string type(void) const;
  56 + std::string name(void) const;
  57 + double solidangle(void) const;
  58 +};
  59 +
  60 +
  61 +/***********************************************************************//**
  62 + * @brief GSkyRegion class extension
  63 + ***************************************************************************/
... ...
pyext/GSkyRegionCircle.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegionCircle.i - Sky region class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GSkyRegionCircle.i
  23 + * @brief Sky region class SWIG file.
  24 + * @author Michael Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GSkyRegionCircle.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GSkyRegionCircle
  35 + *
  36 + * @brief GSkyRegionCircle dervied sky region class
  37 + ***************************************************************************/
  38 +class GSkyRegionCircle : public GSkyRegion {
  39 +public:
  40 + // Constructors and destructors
  41 + GSkyRegionCircle(void);
  42 + GSkyRegionCircle(const GSkyRegionCircle& circle);
  43 + explicit GSkyRegionCircle(const std::string &line);
  44 + GSkyRegionCircle(GSkyDir &centre, const double &radius);
  45 + // GSkyRegionCircle(const double &ra, const double &dec, const double &radius);
  46 + virtual ~GSkyRegionCircle(void);
  47 +
  48 + // Implemented methods
  49 + void clear(void);
  50 + GSkyRegionCircle* clone(void) const;
  51 + GSkyDir centre(void) const;
  52 + double radius(void) const;
  53 + void radius(const double& radius);
  54 + void centre(const GSkyDir& centre);
  55 + void centre(const double& ra,const double& dec);
  56 + void read(const std::string line);
  57 + std::string write(void) const;
  58 + bool contains(const GSkyDir& dir) const;
  59 + bool contains(const GSkyRegion& reg) const;
  60 + bool overlaps(const GSkyRegion& reg) const;
  61 +};
  62 +
  63 +
  64 +/***********************************************************************//**
  65 + * @brief GSkyRegionCircle class extension
  66 + ***************************************************************************/
  67 +%extend GSkyRegionCircle {
  68 + GSkyRegionCircle copy() {
  69 + return (*self);
  70 + }
  71 +};
... ...
pyext/GSkyRegions.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegions.i - Sky region container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GSkyRegions.i
  23 + * @brief Sky regions container class SWIG file
  24 + * @author Pierrick Martin
  25 + */
  26 +
  27 +%{
  28 +/* Put headers and other declarations here that are needed for compilation */
  29 +#include "GBase.hpp"
  30 +#include "GSkyRegion.hpp"
  31 +%}
  32 +
  33 +/***********************************************************************//**
  34 + * @class GSkyRegions
  35 + *
  36 + * @brief Sky region container class
  37 + ***************************************************************************/
  38 +class GSkyRegions : public GBase {
  39 +public:
  40 + // Constructors and destructors
  41 + GSkyRegions(void);
  42 + GSkyRegions(const GSkyRegions& regions);
  43 + explicit GSkyRegions(const std::string& filename);
  44 + virtual ~GSkyRegions(void);
  45 +
  46 + // Methods
  47 + void clear(void);
  48 + GSkyRegions* clone(void) const;
  49 + GSkyRegion* at(const int& index);
  50 + const GSkyRegion* at(const int& index) const;
  51 + int size(void) const;
  52 + bool isempty(void) const;
  53 + GSkyRegion* set(const int& index, const GSkyRegion& region);
  54 + GSkyRegion* set(const std::string& name, const GSkyRegion& region);
  55 + GSkyRegion* append(const GSkyRegion& region);
  56 + GSkyRegion* insert(const int& index, const GSkyRegion& region);
  57 + GSkyRegion* insert(const std::string& name, const GSkyRegion& region);
  58 + void remove(const int& index);
  59 + void remove(const std::string& name);
  60 + void reserve(const int& num);
  61 + void extend(const GSkyRegions& regions);
  62 + bool contains(const std::string& name) const;
  63 + void load(const std::string& filename);
  64 + void save(const std::string& filename) const;
  65 +
  66 +};
  67 +
  68 +/***********************************************************************//**
  69 + * @brief GSkyRegions class extension
  70 + ***************************************************************************/
  71 +%extend GSkyRegions {
  72 + GSkyRegion* __getitem__(const int& index) {
  73 + if (index >= 0 && index < self->size()) {
  74 + return (*self)[index];
  75 + }
  76 + else {
  77 + throw GException::out_of_range("__getitem__(int)", index, self->size());
  78 + }
  79 + }
  80 + void __setitem__(const int& index, const GSkyRegion& val) {
  81 + if (index>=0 && index < self->size()) {
  82 + self->set(index, val);
  83 + return;
  84 + }
  85 + else {
  86 + throw GException::out_of_range("__setitem__(int)", index, self->size());
  87 + }
  88 + }
  89 + int __len__() {
  90 + return (self->size());
  91 + }
  92 + GSkyRegions copy() {
  93 + return (*self);
  94 + }
  95 +};
... ...
pyext/Makefile.am
  View file @7656e43
... ... @@ -62,6 +62,7 @@ BUILT_SOURCES = gammalib/app_wrap.cpp gammalib/app.py \
62 62 gammalib/sky_wrap.cpp gammalib/sky.py \
63 63 gammalib/support_wrap.cpp gammalib/support.py \
64 64 gammalib/xml_wrap.cpp gammalib/xml.py \
  65 + gammalib/xspec_wrap.cpp gammalib/xspec.py \
65 66 gammalib/test_wrap.cpp gammalib/test.py \
66 67 $(INST_SWIG_TARGETS)
67 68  
... ...
pyext/gammalib/sky.i
  View file @7656e43
... ... @@ -50,6 +50,8 @@
50 50 %include "GSkyDir.i"
51 51 %include "GSkyPixel.i"
52 52 %include "GSkymap.i"
  53 +%include "GSkyRegion.i"
  54 +%include "GSkyRegionCircle.i"
53 55 %include "GWcs.i"
54 56 %include "GWcsRegistry.i"
55 57 %include "GWcslib.i"
... ...
pyext/gammalib/xspec.i 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * xspec.i - XSPEC module *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + * ----------------------------------------------------------------------- *
  21 + * Usage: *
  22 + * swig -c++ -python -Wall xspec.i *
  23 + ***************************************************************************/
  24 +/**
  25 + * @file xspec.i
  26 + * @brief XSPEC module
  27 + * @author Juergen Knoedlseder
  28 + */
  29 +%module xspec
  30 +%feature("autodoc", "1");
  31 +
  32 +/* __ Headers needed for compilation _____________________________________ */
  33 +%{
  34 +#include <stddef.h>
  35 +#include "GException.hpp"
  36 +#include "GTools.hpp"
  37 +%}
  38 +
  39 +/* __ Include standard typemaps for vectors and strings __________________ */
  40 +%include stl.i
  41 +
  42 +/* __ Include interface classes __________________________________________ */
  43 +%import(module="gammalib.base") "GBase.i";
  44 +//%import(module="gammalib.base") "GContainer.i";
  45 +
  46 +/* __ Make sure that exceptions are catched ______________________________ */
  47 +%import(module="gammalib.support") "GException.i";
  48 +
  49 +/* __ Application module _________________________________________________ */
  50 +%include "GArf.i"
  51 +%include "GPha.i"
  52 +%include "GRmf.i"
... ...
pyext/setup.py.in
  View file @7656e43
1 1 #############################################################################
2 2 # Disutils setup file for gammalib Python extension module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2011-2012 Juergen Knoedlseder #
  4 +# Copyright (C) 2011-2013 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -31,7 +31,7 @@ pkg_name = &#39;gammalib&#39;
31 31  
32 32 # Define core modules
33 33 core_modules = ['app', 'base', 'fits', 'linalg', 'model', 'numerics', \
34   - 'obs', 'opt', 'sky', 'support', 'test', 'xml']
  34 + 'obs', 'opt', 'sky', 'support', 'test', 'xml', 'xspec']
35 35  
36 36 # Detect execution mode
37 37 is_build = True
... ...
src/Makefile.am
  View file @7656e43
... ... @@ -16,6 +16,7 @@ SUBDIRS = base \
16 16 fits \
17 17 xml \
18 18 vo \
  19 + xspec \
19 20 sky \
20 21 opt \
21 22 obs \
... ... @@ -59,6 +60,7 @@ libgamma_la_LIBADD = base/libbase.la \
59 60 fits/libfits.la \
60 61 xml/libxml.la \
61 62 vo/libvo.la \
  63 + xspec/libxspec.la \
62 64 sky/libsky.la \
63 65 opt/libopt.la \
64 66 obs/libobs.la \
... ...
src/obs/GEbounds.cpp
  View file @7656e43
... ... @@ -699,8 +699,14 @@ void GEbounds::write(GFits* file, const std::string&amp; extname) const
699 699 * @param[in] eng Energy.
700 700 * @return Bin index.
701 701 *
702   - * Returns the energy boundary bin index for a given energy. If the energy
703   - * falls outside all boundaries, -1 is returned.
  702 + * Returns the energy boundary bin index for a given energy. By convention,
  703 + * the limits for an energy bin are defined as
  704 + *
  705 + * min <= energy < max
  706 + *
  707 + * i.e. and energy equals to max falls above the largest energy.
  708 + *
  709 + * If the energy falls outside all boundaries, -1 is returned.
704 710 ***************************************************************************/
705 711 int GEbounds::index(const GEnergy& eng) const
706 712 {
... ... @@ -709,7 +715,7 @@ int GEbounds::index(const GEnergy&amp; eng) const
709 715  
710 716 // Search all energy boundaries for containment
711 717 for (int i = 0; i < m_num; ++i) {
712   - if (eng >= m_min[i] && eng <= m_max[i]) {
  718 + if (eng >= m_min[i] && eng < m_max[i]) {
713 719 index = i;
714 720 break;
715 721 }
... ...
src/sky/GSkyRegion.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegion.cpp - Abstract virtual sky region base class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GSkyRegion.cpp
  23 + * @brief Abstract sky region base class interface definition
  24 + * @author Pierrick Martin
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GSkyRegion.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +
  37 +/* __ Macros _____________________________________________________________ */
  38 +
  39 +/* __ Coding definitions _________________________________________________ */
  40 +
  41 +/* __ Debug definitions __________________________________________________ */
  42 +
  43 +
  44 +/*==========================================================================
  45 + = =
  46 + = Constructors/destructors =
  47 + = =
  48 + ==========================================================================*/
  49 +
  50 +/***********************************************************************//**
  51 + * @brief Void constructor
  52 + ***************************************************************************/
  53 +GSkyRegion::GSkyRegion(void)
  54 +{
  55 + // Initialise members
  56 + init_members();
  57 +
  58 + // Return
  59 + return;
  60 +}
  61 +
  62 +
  63 +/***********************************************************************//**
  64 + * @brief Copy constructor
  65 + *
  66 + * @param[in] region sky region.
  67 + ***************************************************************************/
  68 +GSkyRegion::GSkyRegion(const GSkyRegion& region)
  69 +{
  70 + // Initialise members
  71 + init_members();
  72 +
  73 + // Copy members
  74 + copy_members(region);
  75 +
  76 + // Return
  77 + return;
  78 +}
  79 +
  80 +
  81 +/***********************************************************************//**
  82 + * @brief Destructor
  83 + ***************************************************************************/
  84 +GSkyRegion::~GSkyRegion(void)
  85 +{
  86 + // Free members
  87 + free_members();
  88 +
  89 + // Return
  90 + return;
  91 +}
  92 +
  93 +
  94 +/*==========================================================================
  95 + = =
  96 + = Operators =
  97 + = =
  98 + ==========================================================================*/
  99 +
  100 +/***********************************************************************//**
  101 + * @brief Assignment operator
  102 + *
  103 + * @param[in] region sky region.
  104 + * @return sky region.
  105 + ***************************************************************************/
  106 +GSkyRegion& GSkyRegion::operator=(const GSkyRegion& region)
  107 +{
  108 + // Execute only if object is not identical
  109 + if (this != &region) {
  110 +
  111 + // Free members
  112 + free_members();
  113 +
  114 + // Initialise private members for clean destruction
  115 + init_members();
  116 +
  117 + // Copy members
  118 + copy_members(region);
  119 +
  120 + } // endif: object was not identical
  121 +
  122 + // Return
  123 + return *this;
  124 +}
  125 +
  126 +
  127 +/*==========================================================================
  128 + = =
  129 + = Private methods =
  130 + = =
  131 + ==========================================================================*/
  132 +
  133 +/***********************************************************************//**
  134 + * @brief Initialise class members
  135 + ***************************************************************************/
  136 +void GSkyRegion::init_members(void)
  137 +{
  138 + // Initialise members
  139 + m_type.clear();
  140 + m_name.clear();
  141 + m_solid = 0.0;
  142 +
  143 + // Return
  144 + return;
  145 +}
  146 +
  147 +
  148 +/***********************************************************************//**
  149 + * @brief Copy class members
  150 + *
  151 + * @param[in] region sky region.
  152 + ***************************************************************************/
  153 +void GSkyRegion::copy_members(const GSkyRegion& region)
  154 +{
  155 + // Copy members
  156 + m_type = region.m_type;
  157 + m_name = region.m_name;
  158 + m_solid = region.m_solid;
  159 +
  160 + // Return
  161 + return;
  162 +}
  163 +
  164 +
  165 +/***********************************************************************//**
  166 + * @brief Delete class members
  167 + ***************************************************************************/
  168 +void GSkyRegion::free_members(void)
  169 +{
  170 + // Return
  171 + return;
  172 +}
... ...
src/sky/GSkyRegionCircle.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegionCircle.cpp - Circular sky region class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GSkyRegionCircle.hpp
  23 + * @brief Circular sky region implementation
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GSkyRegionCircle.hpp"
  32 +#include "GTools.hpp"
  33 +
  34 +/* __ Method name definitions ____________________________________________ */
  35 +#define G_CONSTRUCTOR "GSkyRegionCircle::GSkyRegionCircle(GSkyDir&, double&)"
  36 +#define G_RADIUS "GSkyRegionCircle::radius(double&)"
  37 +#define G_READ "GSkyRegionCircle::read(std::string&)"
  38 +#define G_CONTAINS "GSkyRegionCircle::contains(GSkyRegion&)"
  39 +#define G_OVERLAPS "GSkyRegionCircle::overlaps(GSkyRegion&)"
  40 +
  41 +/* __ Macros _____________________________________________________________ */
  42 +
  43 +/* __ Coding definitions _________________________________________________ */
  44 +
  45 +/* __ Debug definitions __________________________________________________ */
  46 +
  47 +/* __ Prototype __________________________________________________________ */
  48 +
  49 +
  50 +/*==========================================================================
  51 + = =
  52 + = Constructors/destructors =
  53 + = =
  54 + ==========================================================================*/
  55 +
  56 +/***********************************************************************//**
  57 + * @brief Void constructor
  58 + ***************************************************************************/
  59 +GSkyRegionCircle::GSkyRegionCircle(void)
  60 +{
  61 + // Initialise members
  62 + init_members();
  63 +
  64 + // Return
  65 + return;
  66 +}
  67 +
  68 +
  69 +/***********************************************************************//**
  70 + * @brief Direction constructor
  71 + *
  72 + * @param[in] centre GSkyDir of Circle center.
  73 + * @param[in] radius Radius of the region [deg].
  74 + *
  75 + * @exception GException::invalid_argument
  76 + * Radius value is less than 0.
  77 + ***************************************************************************/
  78 +GSkyRegionCircle::GSkyRegionCircle(GSkyDir& centre, const double& radius)
  79 +{
  80 + // Initialise members
  81 + init_members();
  82 +
  83 + // set members
  84 + m_centre = centre;
  85 + m_radius = radius;
  86 +
  87 + // Check if radius is valid
  88 + if (m_radius < 0.0) {
  89 + throw GException::invalid_argument(G_CONSTRUCTOR,
  90 + "Radius of a region can't be less than 0");
  91 + }
  92 +
  93 + // compute solid angle
  94 + compute_solid();
  95 +
  96 + // Return
  97 + return;
  98 +}
  99 +
  100 +/***********************************************************************//**
  101 + * @brief string constructor
  102 + *
  103 + * @param[in] line a string containing a line of a ds9 region file
  104 + ***************************************************************************/
  105 +GSkyRegionCircle::GSkyRegionCircle(const std::string& line)
  106 +{
  107 + // Initialise members
  108 + init_members();
  109 +
  110 + // read from line
  111 + read(line);
  112 +
  113 + // Return
  114 + return;
  115 +
  116 +}
  117 +
  118 +
  119 +/***********************************************************************//**
  120 + * @brief Copy constructor
  121 + *
  122 + * @param[in] circle Circular sky region.
  123 + ***************************************************************************/
  124 +GSkyRegionCircle::GSkyRegionCircle(const GSkyRegionCircle& circle)
  125 +{
  126 + // Initialise members
  127 + init_members();
  128 +
  129 + // Copy members
  130 + copy_members(circle);
  131 +
  132 + // Return
  133 + return;
  134 +}
  135 +
  136 +
  137 +/***********************************************************************//**
  138 + * @brief Destructor
  139 + ***************************************************************************/
  140 +GSkyRegionCircle::~GSkyRegionCircle(void)
  141 +{
  142 + // Free members
  143 + free_members();
  144 +
  145 + // Return
  146 + return;
  147 +}
  148 +
  149 +
  150 +/*==========================================================================
  151 + = =
  152 + = Operators =
  153 + = =
  154 + ==========================================================================*/
  155 +
  156 +/***********************************************************************//**
  157 + * @brief Assignment operator
  158 + *
  159 + * @param[in] circle Circular sky region.
  160 + * @return Circular sky region.
  161 + ***************************************************************************/
  162 +GSkyRegionCircle& GSkyRegionCircle::operator= (const GSkyRegionCircle& circle)
  163 +{
  164 + // Execute only if object is not identical
  165 + if (this != &circle) {
  166 +
  167 + // Free members
  168 + free_members();
  169 +
  170 + // Initialise private members for clean destruction
  171 + init_members();
  172 +
  173 + // Copy members
  174 + copy_members(circle);
  175 +
  176 + } // endif: object was not identical
  177 +
  178 + // Return this object
  179 + return *this;
  180 +}
  181 +
  182 +
  183 +/*==========================================================================
  184 + = =
  185 + = Public methods =
  186 + = =
  187 + ==========================================================================*/
  188 +
  189 +/***********************************************************************//**
  190 + * @brief Clear instance
  191 + ***************************************************************************/
  192 +void GSkyRegionCircle::clear(void)
  193 +{
  194 + // Free members
  195 + free_members();
  196 + this->GSkyRegion::free_members();
  197 +
  198 + // Initialise private members
  199 + this->GSkyRegion::init_members();
  200 + init_members();
  201 +
  202 + // Return
  203 + return;
  204 +}
  205 +
  206 +
  207 +/***********************************************************************//**
  208 + * @brief Clone object
  209 + ***************************************************************************/
  210 +GSkyRegionCircle* GSkyRegionCircle::clone(void) const
  211 +{
  212 + // Clone this Region
  213 + return new GSkyRegionCircle(*this);
  214 +}
  215 +
  216 +
  217 +/***********************************************************************//**
  218 + * @brief Set radius value of region
  219 + *
  220 + * @param[in] radius Radius value [deg].
  221 + *
  222 + * @exception GException::invalid_argument
  223 + * Radius value is less than 0.
  224 + ***************************************************************************/
  225 +void GSkyRegionCircle::radius(const double& radius)
  226 +{
  227 + // Set radius value
  228 + m_radius = radius;
  229 +
  230 + if (m_radius < 0.0) {
  231 + throw GException::invalid_argument(G_RADIUS,
  232 + "Radius of a region can't be less than 0");
  233 + }
  234 +
  235 + // Recompute the solid angle
  236 + compute_solid();
  237 +}
  238 +
  239 +
  240 +/***********************************************************************//**
  241 + * @brief Read region from string
  242 + *
  243 + * @param[in] line String in ds9 format.
  244 + *
  245 + * @exception GException::invalid_value
  246 + * Invalid value found in ds9 format string.
  247 + ***************************************************************************/
  248 +void GSkyRegionCircle::read(const std::string& line)
  249 +{
  250 + // Clear the current instance
  251 + clear();
  252 +
  253 + // Split the string into 2 parts seperated by #
  254 + std::vector<std::string> substrings = gammalib::split(line,"#");
  255 + std::string region_def = (substrings.size() > 0) ? substrings[0] : "";
  256 + std::string comment = (substrings.size() > 1) ? substrings[1] : "";
  257 +
  258 + // Finding the circle
  259 + if (region_def.find("circle") == std::string::npos) {
  260 + throw GException::invalid_value(G_READ,
  261 + "Cannnot find the key word \"circle\" in provided string");
  262 + }
  263 +
  264 + // Get the the coordinate system of the values
  265 + std::string system = gammalib::split(region_def,";")[0];
  266 +
  267 + // Get the substring of the important values
  268 + unsigned pos = region_def.find("circle(");
  269 + unsigned end = region_def.find(")");
  270 + std::string circlestring = region_def.substr(pos+7,end);
  271 + circlestring.erase(circlestring.find(")"),1);
  272 +
  273 + // Get the values of the region x,y,and radius
  274 + std::vector<std::string> values = gammalib::split(circlestring,",");
  275 + if (values.size()!=3) {
  276 + throw GException::invalid_value(G_READ,
  277 + "Invalid number of arguments, circle takes exactly 3 arguments");
  278 + }
  279 + double x = gammalib::todouble(values[0]);
  280 + double y = gammalib::todouble(values[1]);
  281 + double radius = gammalib::todouble(values[2]);
  282 +
  283 + // Initialise sky dir
  284 + GSkyDir dir = GSkyDir();
  285 +
  286 + // Set the values in the correct system
  287 + if (system == "fk5") {
  288 + dir.radec_deg(x,y);
  289 + }
  290 + else if (system == "galactic") {
  291 + dir.lb_deg(x,y);
  292 + }
  293 + else {
  294 + std::stringstream s;
  295 + s << "Provided coordinate system \"" << system << "\" is not supported.";
  296 + throw GException::invalid_value(G_READ, s.str());
  297 + }
  298 +
  299 + // Setting par values
  300 + m_centre = dir;
  301 + m_radius = radius;
  302 +
  303 + // Throwing exception if radius is less than 0
  304 + if (m_radius < 0.0) {
  305 + throw GException::invalid_value(G_READ,
  306 + "Radius of a region can't be less than 0");
  307 + }
  308 +
  309 + // Compute solid angle
  310 + compute_solid();
  311 +
  312 + // Check if there is a given name for the region and set it
  313 + std::vector<std::string>comments = gammalib::split(comment, " ");
  314 + for (int i = 0; i < comments.size(); i++) {
  315 + if (gammalib::contains(comments[i], "text")) {
  316 + std::vector<std::string> attributes = gammalib::split(comments[i], "=");
  317 + if (attributes.size() < 2) {
  318 + throw GException::invalid_value(G_READ,
  319 + "Invalid number of arguments, type of attribute must"
  320 + " be key=value");
  321 + }
  322 + m_name = attributes[1];
  323 + }
  324 + }
  325 +
  326 + // Return
  327 + return;
  328 +}
  329 +
  330 +
  331 +/***********************************************************************//**
  332 + * @brief write region into a string
  333 + *
  334 + * @return String to be written in a ds9 region file
  335 + ***************************************************************************/
  336 +std::string GSkyRegionCircle::write(void) const
  337 +{
  338 + std::string result;
  339 + result.append("fk5;circle(");
  340 + result.append(gammalib::str(m_centre.ra_deg()));
  341 + result.append(",");
  342 + result.append(gammalib::str(m_centre.dec_deg()));
  343 + result.append(",");
  344 + result.append(gammalib::str(m_radius));
  345 + result.append(") # text=");
  346 + result.append(m_name);
  347 + result.append("\n");
  348 + return result;
  349 +}
  350 +
  351 +
  352 +/***********************************************************************//**
  353 + * @brief Print circular region
  354 + *
  355 + * @param[in] chatter Chattiness (defaults to NORMAL).
  356 + * @return String containing region information.
  357 + ***************************************************************************/
  358 +std::string GSkyRegionCircle::print(const GChatter& chatter) const
  359 +{
  360 + // Initialise result string
  361 + std::string result;;
  362 +
  363 + // Continue only if chatter is not silent
  364 + if (chatter != SILENT) {
  365 +
  366 + // Append string
  367 + result.append("=== GSkyRegionCircle ===");
  368 + result.append("\n(");
  369 + result.append(gammalib::str(m_centre.ra_deg()));
  370 + result.append(",");
  371 + result.append(gammalib::str(m_centre.dec_deg()));
  372 + result.append(",");
  373 + result.append(gammalib::str(m_radius));
  374 + result.append(")");
  375 +
  376 + } // endif: chatter was not silent
  377 +
  378 + // Return result
  379 + return result;
  380 +}
  381 +
  382 +/***********************************************************************//**
  383 + * @brief checks if coordinate is within region
  384 + *
  385 + * @param[in] dir Sky direction.
  386 + ***************************************************************************/
  387 +bool GSkyRegionCircle::contains(const GSkyDir& dir) const
  388 +{
  389 + // Initialise return value
  390 + bool dir_is_in = false;
  391 +
  392 + // calculate distance from dir to centre
  393 + double distance = dir.dist_deg(m_centre);
  394 +
  395 + // Check if distance < radius
  396 + if (distance <= radius()) {
  397 + dir_is_in = true;
  398 + }
  399 +
  400 + // Return bool
  401 + return dir_is_in;
  402 +}
  403 +
  404 +/***********************************************************************//**
  405 + * @brief checks if region is fully contained within this region
  406 + *
  407 + * @param[in] reg Sky region.
  408 + *
  409 + * @exception GException::feature_not_implemented
  410 + * Regions differ in type.
  411 + ***************************************************************************/
  412 +bool GSkyRegionCircle::contains(const GSkyRegion& reg) const
  413 +{
  414 + // Initialise return value
  415 + bool fully_inside = false;
  416 +
  417 + // If other region is circle use a simple way to calculate
  418 + if (reg.type() == "Circle") {
  419 +
  420 + // create circular region from reg
  421 + GSkyRegionCircle* regcirc = (GSkyRegionCircle*)reg.clone();
  422 +
  423 + // calculate angular distance between the centers
  424 + double ang_dist = m_centre.dist_deg(regcirc->centre());
  425 +
  426 + // Check if the region is contained in this
  427 + if ((ang_dist + regcirc->radius()) <= m_radius) {
  428 + fully_inside = true;
  429 + }
  430 +
  431 + }
  432 +
  433 + // ... otherwise throw an exception
  434 + else {
  435 + throw GException::feature_not_implemented(G_CONTAINS,
  436 + "Cannot compare two different region types yet");
  437 + }
  438 +
  439 + // Return value
  440 + return fully_inside;
  441 +}
  442 +
  443 +/***********************************************************************//**
  444 + * @brief Checks if region is overlapping with this region
  445 + *
  446 + * @param[in] reg Sky region.
  447 + *
  448 + * @exception GException::feature_not_implemented
  449 + * Regions differ in type.
  450 + ***************************************************************************/
  451 +bool GSkyRegionCircle::overlaps(const GSkyRegion& reg) const
  452 +{
  453 + // Initialise return value
  454 + bool overlap = false;
  455 +
  456 + // If other region is circle use a simple way to calculate
  457 + if (reg.type() == "Circle") {
  458 +
  459 + // Create circular region from reg
  460 + GSkyRegionCircle* regcirc = (GSkyRegionCircle*)reg.clone();
  461 +
  462 + // Calculate angular distance between the centers
  463 + double ang_dist = m_centre.dist_deg(regcirc->centre());
  464 +
  465 + // Check if the distance is smaller than the sum of both radii
  466 + if (ang_dist <= (m_radius + regcirc->radius())) {
  467 + overlap = true;
  468 + }
  469 +
  470 + }
  471 +
  472 + // ... otherwise throw an exception
  473 + else {
  474 + throw GException::feature_not_implemented(G_OVERLAPS,
  475 + "Cannot compare two different region types yet");
  476 + }
  477 +
  478 + // Return value
  479 + return overlap;
  480 +}
  481 +
  482 +
  483 +
  484 +/*==========================================================================
  485 + = =
  486 + = Private methods =
  487 + = =
  488 + ==========================================================================*/
  489 +
  490 +/***********************************************************************//**
  491 + * @brief Initialise class members
  492 + ***************************************************************************/
  493 +void GSkyRegionCircle::init_members(void)
  494 +{
  495 + // Initialise members
  496 + m_centre = GSkyDir();
  497 + m_radius = 0.0;
  498 + m_type = "Circle";
  499 +
  500 + //Return
  501 + return;
  502 +}
  503 +
  504 +
  505 +/***********************************************************************//**
  506 + * @brief Copy class members
  507 + *
  508 + * @param[in] region Circular sky region.
  509 + ***************************************************************************/
  510 +void GSkyRegionCircle::copy_members(const GSkyRegionCircle& region)
  511 +{
  512 + // Copy attributes
  513 + m_centre = region.m_centre;
  514 + m_radius = region.m_radius;
  515 +
  516 + // Return
  517 + return;
  518 +}
  519 +
  520 +
  521 +/***********************************************************************//**
  522 + * @brief Delete class members
  523 + ***************************************************************************/
  524 +void GSkyRegionCircle::free_members(void)
  525 +{
  526 + // Return
  527 + return;
  528 +}
  529 +
  530 +
  531 +/***********************************************************************//**
  532 + * @brief compute solid angle
  533 + ***************************************************************************/
  534 +void GSkyRegionCircle::compute_solid(void)
  535 +{
  536 + // convert radius to radians
  537 + double radius_rad = m_radius * gammalib::deg2rad;
  538 +
  539 + // compute solid angle
  540 + m_solid = gammalib::twopi * (1.0 - std::cos(radius_rad));
  541 +
  542 + // Return
  543 + return;
  544 +}
... ...
src/sky/GSkyRegions.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GSkyRegions.cpp - Sky region container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 +* @file GSkyRegions.cpp
  23 +* @brief Sky regions container class definition
  24 +* @author Pierrick Martin
  25 +*/
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#include <fstream>
  31 +#endif
  32 +#include "GBase.hpp"
  33 +#include "GSkyRegion.hpp"
  34 +#include "GSkyRegionCircle.hpp"
  35 +#include "GSkyRegions.hpp"
  36 +#include "GTools.hpp"
  37 +
  38 +/* __ Method name definitions ____________________________________________ */
  39 +#define G_ACCESS "GSkyRegions::operator[](std::string&)"
  40 +#define G_AT "GSkyRegions::at(int&)"
  41 +#define G_SET1 "GSkyRegions::set(int&, GSkyRegion&)"
  42 +#define G_SET2 "GSkyRegions::set(std::string&, GSkyRegion&)"
  43 +#define G_APPEND "GSkyRegions::append(GSkyRegion&)"
  44 +#define G_INSERT1 "GSkyRegions::insert(int&, GSkyRegion&)"
  45 +#define G_INSERT2 "GSkyRegions::insert(std::string&, GSkyRegion&)"
  46 +#define G_REMOVE1 "GSkyRegions::remove(int&)"
  47 +#define G_REMOVE2 "GSkyRegions::remove(std::string&)"
  48 +#define G_EXTEND "GSkyRegions::extend(GSkyRegions&)"
  49 +#define G_LOAD "GSkyRegions::load(std::string&)"
  50 +#define G_SAVE "GSkyRegions::save(std::string&)"
  51 +
  52 +/* __ Macros _____________________________________________________________ */
  53 +
  54 +/* __ Coding definitions _________________________________________________ */
  55 +
  56 +/* __ Debug definitions __________________________________________________ */
  57 +
  58 +
  59 +/*==========================================================================
  60 + = =
  61 + = Constructors/destructors =
  62 + = =
  63 + ==========================================================================*/
  64 +
  65 +/***********************************************************************//**
  66 + * @brief Void constructor
  67 + ***************************************************************************/
  68 +GSkyRegions::GSkyRegions(void)
  69 +{
  70 + // Initialise members
  71 + init_members();
  72 +
  73 + // Return
  74 + return;
  75 +}
  76 +
  77 +
  78 +/***********************************************************************//**
  79 + * @brief Copy constructor
  80 + *
  81 + * @param[in] regions region container.
  82 + ***************************************************************************/
  83 +GSkyRegions::GSkyRegions(const GSkyRegions& regions)
  84 +{
  85 + // Initialise members
  86 + init_members();
  87 +
  88 + // Copy members
  89 + copy_members(regions);
  90 +
  91 + // Return
  92 + return;
  93 +}
  94 +
  95 +
  96 +/***********************************************************************//**
  97 + * @brief Load constructor
  98 + *
  99 + * @param[in] filename DS9 region file.
  100 + *
  101 + * Constructs region container from a DS9 region file.
  102 + ***************************************************************************/
  103 +GSkyRegions::GSkyRegions(const std::string& filename)
  104 +{
  105 + // Initialise members
  106 + init_members();
  107 +
  108 + // Load XML file
  109 + load(filename);
  110 +
  111 + // Return
  112 + return;
  113 +}
  114 +
  115 +
  116 +/***********************************************************************//**
  117 + * @brief Destructor
  118 + ***************************************************************************/
  119 +GSkyRegions::~GSkyRegions(void)
  120 +{
  121 + // Free members
  122 + free_members();
  123 +
  124 + // Return
  125 + return;
  126 +}
  127 +
  128 +
  129 +/*==========================================================================
  130 + = =
  131 + = Operators =
  132 + = =
  133 + ==========================================================================*/
  134 +
  135 +/***********************************************************************//**
  136 + * @brief Assignment operator
  137 + *
  138 + * @param[in] regions region container.
  139 + * @return region container.
  140 + ***************************************************************************/
  141 +GSkyRegions& GSkyRegions::operator=(const GSkyRegions& regions)
  142 +{
  143 + // Execute only if object is not identical
  144 + if (this != &regions) {
  145 +
  146 + // Free members
  147 + free_members();
  148 +
  149 + // Initialise members
  150 + init_members();
  151 +
  152 + // Copy members
  153 + copy_members(regions);
  154 +
  155 + } // endif: object was not identical
  156 +
  157 + // Return
  158 + return *this;
  159 +}
  160 +
  161 +
  162 +/***********************************************************************//**
  163 + * @brief Return pointer to region
  164 + *
  165 + * @param[in] name region name.
  166 + *
  167 + * @exception GException::invalid_argument
  168 + * region with specified name not found in container.
  169 + *
  170 + * Returns a pointer to the region with the specified @p name.
  171 + ***************************************************************************/
  172 +GSkyRegion* GSkyRegions::operator[](const std::string& name)
  173 +{
  174 + // Get region index
  175 + int index = get_index(name);
  176 +
  177 + // Throw exception if region name was not found
  178 + if (index == -1) {
  179 + std::string msg ="Region with name \""+name+"\" not found in container.";
  180 + throw GException::invalid_argument(G_ACCESS, msg);
  181 + }
  182 +
  183 + // Return pointer
  184 + return m_regions[index];
  185 +}
  186 +
  187 +
  188 +/***********************************************************************//**
  189 + * @brief Return pointer to region (const version)
  190 + *
  191 + * @param[in] name region name.
  192 + *
  193 + * @exception GException::invalid_argument
  194 + * region with specified name not found in container.
  195 + *
  196 + * Returns a const pointer to the region with the specified @p name.
  197 + ***************************************************************************/
  198 +const GSkyRegion* GSkyRegions::operator[](const std::string& name) const
  199 +{
  200 + // Get region index
  201 + int index = get_index(name);
  202 +
  203 + // Throw exception if region name was not found
  204 + if (index == -1) {
  205 + std::string msg ="Region with name \""+name+"\" not found in container.";
  206 + throw GException::invalid_argument(G_ACCESS, msg);
  207 + }
  208 +
  209 + // Return pointer
  210 + return m_regions[index];
  211 +}
  212 +
  213 +
  214 +/*==========================================================================
  215 + = =
  216 + = Public methods =
  217 + = =
  218 + ==========================================================================*/
  219 +
  220 +/***********************************************************************//**
  221 + * @brief Clear object
  222 + *
  223 + * Removes all regions from the container.
  224 + ***************************************************************************/
  225 +void GSkyRegions::clear(void)
  226 +{
  227 + // Free class members
  228 + free_members();
  229 +
  230 + // Initialise members
  231 + init_members();
  232 +
  233 + // Return
  234 + return;
  235 +}
  236 +
  237 +
  238 +/***********************************************************************//**
  239 + * @brief Clone instance
  240 + *
  241 + * @return Pointer to deep copy of region container
  242 + *
  243 + * Makes a deep copy of the region container instance.
  244 + ***************************************************************************/
  245 +GSkyRegions* GSkyRegions::clone(void) const
  246 +{
  247 + return new GSkyRegions(*this);
  248 +}
  249 +
  250 +
  251 +/***********************************************************************//**
  252 + * @brief Return pointer to region
  253 + *
  254 + * @param[in] index region index [0,...,size()-1].
  255 + *
  256 + * @exception GException::out_of_range
  257 + * region index is out of range.
  258 + *
  259 + * Returns a pointer to the region with the specified @p index.
  260 + ***************************************************************************/
  261 +GSkyRegion* GSkyRegions::at(const int& index)
  262 +{
  263 + // Raise exception if index is out of range
  264 + if (index < 0 || index >= size()) {
  265 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  266 + }
  267 +
  268 + // Return pointer
  269 + return m_regions[index];
  270 +}
  271 +
  272 +
  273 +/***********************************************************************//**
  274 + * @brief Return pointer to region (const version)
  275 + *
  276 + * @param[in] index region index [0,...,size()-1].
  277 + *
  278 + * @exception GException::out_of_range
  279 + * region index is out of range.
  280 + *
  281 + * Returns a const pointer to the region with the specified @p index.
  282 + ***************************************************************************/
  283 +const GSkyRegion* GSkyRegions::at(const int& index) const
  284 +{
  285 + // Raise exception if index is out of range
  286 + if (index < 0 || index >= size()) {
  287 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  288 + }
  289 +
  290 + // Return pointer
  291 + return m_regions[index];
  292 +}
  293 +
  294 +
  295 +/***********************************************************************//**
  296 + * @brief Set region in container
  297 + *
  298 + * @param[in] index sky region index [0,...,size()-1].
  299 + * @param[in] region sky region.
  300 + * @return Pointer to deep copy of sky region.
  301 + *
  302 + * @exception GException::out_of_range
  303 + * region index is out of range.
  304 + * @exception GException::invalid_value
  305 + * Name of region exists already in container.
  306 + *
  307 + * Set sky region in the container. A deep copy of the region will be made.
  308 + ***************************************************************************/
  309 +GSkyRegion* GSkyRegions::set(const int& index, const GSkyRegion& region)
  310 +{
  311 + // Compile option: raise exception if index is out of range
  312 + #if defined(G_RANGE_CHECK)
  313 + if (index < 0 || index >= size()) {
  314 + throw GException::out_of_range(G_SET1, index, 0, size()-1);
  315 + }
  316 + #endif
  317 +
  318 + // Check if a region with specified name does not yet exist
  319 + int inx = get_index(region.name());
  320 + if (inx != -1 && inx != index) {
  321 + std::string msg =
  322 + "Attempt to set region with name \""+region.name()+"\" in region"
  323 + " container at index "+gammalib::str(index)+", but a region with"
  324 + " the same name exists already at index "+gammalib::str(inx)+
  325 + " in the container.\n"
  326 + "Every region in the region container needs a unique name.";
  327 + throw GException::invalid_value(G_SET1, msg);
  328 + }
  329 +
  330 + // Delete any existing region
  331 + if (m_regions[index] != NULL) delete m_regions[index];
  332 +
  333 + // Assign new region by cloning
  334 + m_regions[index] = region.clone();
  335 +
  336 + // Return pointer to region
  337 + return m_regions[index];
  338 +}
  339 +
  340 +
  341 +/***********************************************************************//**
  342 + * @brief Set region in container
  343 + *
  344 + * @param[in] name region name.
  345 + * @param[in] region region pointer.
  346 + * @return Pointer to deep copy of region.
  347 + *
  348 + * @exception GException::invalid_argument
  349 + * region with specified name not found in container.
  350 + * @exception GException::invalid_value
  351 + * Name of region exists already in container.
  352 + *
  353 + * Set region in the container. A deep copy of the region will be made.
  354 + ***************************************************************************/
  355 +GSkyRegion* GSkyRegions::set(const std::string& name, const GSkyRegion& region)
  356 +{
  357 + // Get parameter index
  358 + int index = get_index(name);
  359 +
  360 + // Throw exception if parameter name was not found
  361 + if (index == -1) {
  362 + throw GException::invalid_argument(G_SET2, name);
  363 + }
  364 +
  365 + // Check if a region with specified name does not yet exist
  366 + int inx = get_index(region.name());
  367 + if (inx != -1 && inx != index) {
  368 + std::string msg =
  369 + "Attempt to set region with name \""+region.name()+"\" in region"
  370 + " container at index "+gammalib::str(index)+", but a region with"
  371 + " the same name exists already at index "+gammalib::str(inx)+
  372 + " in the container.\n"
  373 + "Every region in the region container needs a unique name.";
  374 + throw GException::invalid_value(G_SET2, msg);
  375 + }
  376 +
  377 + // Delete any existing region
  378 + if (m_regions[index] != NULL) delete m_regions[index];
  379 +
  380 + // Assign new region by cloning
  381 + m_regions[index] = region.clone();
  382 +
  383 + // Return pointer to region
  384 + return m_regions[index];
  385 +}
  386 +
  387 +
  388 +/***********************************************************************//**
  389 + * @brief Append region to container
  390 + *
  391 + * @param[in] region region.
  392 + * @return Pointer to deep copy of region.
  393 + *
  394 + * @exception GException::invalid_value
  395 + * Name of region exists already in container.
  396 + *
  397 + * Appends region to the container by making a deep copy of the region and
  398 + * storing its pointer.
  399 + ***************************************************************************/
  400 +GSkyRegion* GSkyRegions::append(const GSkyRegion& region)
  401 +{
  402 + // Check if a region with specified name does not yet exist
  403 + int inx = get_index(region.name());
  404 + if (inx != -1) {
  405 + std::string msg =
  406 + "Attempt to append region with name \""+region.name()+"\" to region"
  407 + " container, but a region with the same name exists already at"
  408 + " index "+gammalib::str(inx)+" in the container.\n"
  409 + "Every region in the region container needs a unique name.";
  410 + throw GException::invalid_value(G_APPEND, msg);
  411 + }
  412 +
  413 + // Create deep copy of region
  414 + GSkyRegion* ptr = region.clone();
  415 +
  416 + // Append deep copy of region
  417 + m_regions.push_back(ptr);
  418 +
  419 + // Return pointer to region
  420 + return ptr;
  421 +}
  422 +
  423 +
  424 +/***********************************************************************//**
  425 + * @brief Insert region into container
  426 + *
  427 + * @param[in] index region index [0,...,size()-1].
  428 + * @param[in] region region.
  429 + * @return Pointer to deep copy of region.
  430 + *
  431 + * @exception GException::out_of_range
  432 + * region index is out of range.
  433 + * @exception GException::invalid_value
  434 + * Name of region exists already in container.
  435 + *
  436 + * Inserts a @p region into the container before the region with the specified
  437 + * @p index.
  438 + ***************************************************************************/
  439 +GSkyRegion* GSkyRegions::insert(const int& index, const GSkyRegion& region)
  440 +{
  441 + // Compile option: raise exception if index is out of range
  442 + #if defined(G_RANGE_CHECK)
  443 + if (isempty()) {
  444 + if (index > 0) {
  445 + throw GException::out_of_range(G_INSERT1, index, 0, size()-1);
  446 + }
  447 + }
  448 + else {
  449 + if (index < 0 || index >= size()) {
  450 + throw GException::out_of_range(G_INSERT1, index, 0, size()-1);
  451 + }
  452 + }
  453 + #endif
  454 +
  455 + // Check if a region with specified name does not yet exist
  456 + int inx = get_index(region.name());
  457 + if (inx != -1) {
  458 + std::string msg =
  459 + "Attempt to insert region with name \""+region.name()+"\" in region"
  460 + " container before index "+gammalib::str(index)+", but a region"
  461 + " with the same name exists already at index "+gammalib::str(inx)+
  462 + " in the container.\n"
  463 + "Every region in the region container needs a unique name.";
  464 + throw GException::invalid_value(G_INSERT1, msg);
  465 + }
  466 +
  467 + // Create deep copy of region
  468 + GSkyRegion* ptr = region.clone();
  469 +
  470 + // Inserts deep copy of region
  471 + m_regions.insert(m_regions.begin()+index, ptr);
  472 +
  473 + // Return pointer to region
  474 + return ptr;
  475 +}
  476 +
  477 +
  478 +/***********************************************************************//**
  479 + * @brief Insert region into container
  480 + *
  481 + * @param[in] name region name.
  482 + * @param[in] region region.
  483 + * @return Pointer to deep copy of region.
  484 + *
  485 + * @exception GException::invalid_argument
  486 + * region with specified name not found in container.
  487 + * @exception GException::invalid_value
  488 + * Name of region exists already in container.
  489 + *
  490 + * Inserts a @p region into the container before the region with the specified
  491 + * @p name.
  492 + ***************************************************************************/
  493 +GSkyRegion* GSkyRegions::insert(const std::string& name, const GSkyRegion& region)
  494 +{
  495 + // Get parameter index
  496 + int index = get_index(name);
  497 +
  498 + // Throw exception if parameter name was not found
  499 + if (index == -1) {
  500 + throw GException::invalid_argument(G_INSERT2, name);
  501 + }
  502 +
  503 + // Check if a region with specified name does not yet exist
  504 + int inx = get_index(region.name());
  505 + if (inx != -1) {
  506 + std::string msg =
  507 + "Attempt to insert region with name \""+region.name()+"\" in region"
  508 + " container before index "+gammalib::str(index)+", but a region"
  509 + " with the same name exists already at index "+gammalib::str(inx)+
  510 + " in the container.\n"
  511 + "Every region in the region container needs a unique name.";
  512 + throw GException::invalid_value(G_INSERT2, msg);
  513 + }
  514 +
  515 + // Create deep copy of region
  516 + GSkyRegion* ptr = region.clone();
  517 +
  518 + // Inserts deep copy of region
  519 + m_regions.insert(m_regions.begin()+index, ptr);
  520 +
  521 + // Return pointer to region
  522 + return ptr;
  523 +}
  524 +
  525 +
  526 +/***********************************************************************//**
  527 + * @brief Remove region from container
  528 + *
  529 + * @param[in] index region index [0,...,size()-1].
  530 + *
  531 + * @exception GException::out_of_range
  532 + * region index is out of range.
  533 + *
  534 + * Remove region of specified @p index from container.
  535 + ***************************************************************************/
  536 +void GSkyRegions::remove(const int& index)
  537 +{
  538 + // Compile option: raise exception if index is out of range
  539 + #if defined(G_RANGE_CHECK)
  540 + if (index < 0 || index >= size()) {
  541 + throw GException::out_of_range(G_REMOVE1, index, 0, size()-1);
  542 + }
  543 + #endif
  544 +
  545 + // Erase region component from container
  546 + m_regions.erase(m_regions.begin() + index);
  547 +
  548 + // Return
  549 + return;
  550 +}
  551 +
  552 +
  553 +/***********************************************************************//**
  554 + * @brief Remove region from container
  555 + *
  556 + * @param[in] name region name.
  557 + *
  558 + * @exception GException::invalid_argument
  559 + * region with specified name not found in container.
  560 + *
  561 + * Remove region of specified @p name from container.
  562 + ***************************************************************************/
  563 +void GSkyRegions::remove(const std::string& name)
  564 +{
  565 + // Get parameter index
  566 + int index = get_index(name);
  567 +
  568 + // Throw exception if region name was not found
  569 + if (index == -1) {
  570 + throw GException::invalid_argument(G_REMOVE2, name);
  571 + }
  572 +
  573 + // Erase region component from container
  574 + m_regions.erase(m_regions.begin() + index);
  575 +
  576 + // Return
  577 + return;
  578 +}
  579 +
  580 +
  581 +/***********************************************************************//**
  582 + * @brief Append region container
  583 + *
  584 + * @param[in] regions region container.
  585 + *
  586 + * Append region container to the container.
  587 + ***************************************************************************/
  588 +void GSkyRegions::extend(const GSkyRegions& regions)
  589 +{
  590 + // Do nothing if region container is empty
  591 + if (!regions.isempty()) {
  592 +
  593 + // Get size. Note that we extract the size first to avoid an
  594 + // endless loop that arises when a container is appended to
  595 + // itself.
  596 + int num = regions.size();
  597 +
  598 + // Reserve enough space
  599 + reserve(size() + num);
  600 +
  601 + // Loop over all region components and append pointers to deep copies
  602 + for (int i = 0; i < num; ++i) {
  603 +
  604 + // Check if region name does not yet exist
  605 + int inx = get_index(regions[i]->name());
  606 + if (inx != -1) {
  607 + std::string msg =
  608 + "Attempt to append region with name \""+regions[i]->name()+
  609 + "\" to region container, but a region with the same name"
  610 + " exists already at index "+gammalib::str(inx)+" in the"
  611 + " container.\n"
  612 + "Every region in the region container needs a unique name.";
  613 + throw GException::invalid_value(G_EXTEND, msg);
  614 + }
  615 +
  616 + // Append region to container
  617 + m_regions.push_back(regions[i]->clone());
  618 +
  619 + } // endfor: looped over all regions
  620 +
  621 + } // endif: region container was not empty
  622 +
  623 + // Return
  624 + return;
  625 +}
  626 +
  627 +
  628 +/***********************************************************************//**
  629 + * @brief Signals if region name exists
  630 + *
  631 + * @param[in] name region name.
  632 + * @return True if region with specified @p name exists.
  633 + *
  634 + * Searches all region names for a match with the specified @p name. If the
  635 + * specified name has been found, true is returned.
  636 + ***************************************************************************/
  637 +bool GSkyRegions::contains(const std::string& name) const
  638 +{
  639 + // Get region index
  640 + int index = get_index(name);
  641 +
  642 + // Return
  643 + return (index != -1);
  644 +}
  645 +
  646 +
  647 +/***********************************************************************//**
  648 + * @brief Load regions from DS9 region file
  649 + *
  650 + * @param[in] filename DS9 region filename.
  651 + *
  652 + * @exception GException::file_open_error
  653 + * File could not be opened.
  654 + *
  655 + * Loads all regions from a DS9 region file. See the read() method for more
  656 + * information about the expected structure of the XML file.
  657 + ***************************************************************************/
  658 +void GSkyRegions::load(const std::string& filename)
  659 +{
  660 + // Clear any existing regions
  661 + clear();
  662 +
  663 + // Open file. Throw an exception if opening failed.
  664 + std::ifstream ds9file;
  665 + ds9file.open(filename.c_str());
  666 + if (ds9file.is_open()) {
  667 +
  668 + // Loop over file lines
  669 + std::string fileline="";
  670 + std::string coordsys="galactic";
  671 + while ( ds9file.good() ) {
  672 +
  673 + // Read one line
  674 + getline (ds9file,fileline);
  675 +
  676 + // If line is a comment then continue
  677 + if (fileline[0] == '#') {
  678 + continue;
  679 + }
  680 +
  681 + // Check for global definition of coordinate system
  682 + if (std::string::npos != fileline.find("fk5")) {
  683 + coordsys="fk5";
  684 + }
  685 +
  686 + // If region is a circle
  687 + if (std::string::npos != fileline.find("circle")) {
  688 +
  689 + // Create instance of GSkyRegion object
  690 + GSkyRegionCircle region;
  691 +
  692 + // If coordinate system and region defined on the same line
  693 + if ((std::string::npos != fileline.find("fk5")) ||
  694 + (std::string::npos != fileline.find("galactic"))) {
  695 + region.read(fileline);
  696 + append(region);
  697 + }
  698 +
  699 + // else, prepend the coordinate system
  700 + else {
  701 + std::string newfileline=coordsys;
  702 + newfileline.append("; ");
  703 + newfileline.append(fileline);
  704 + region.read(newfileline);
  705 + append(region);
  706 + }
  707 + }
  708 + }
  709 +
  710 + // Close file
  711 + ds9file.close();
  712 +
  713 + // Store filename
  714 + m_filename = filename;
  715 +
  716 + }
  717 +
  718 + // File could not be opened
  719 + else {
  720 + throw GException::file_open_error(G_LOAD, filename);
  721 + }
  722 +
  723 + // Return
  724 + return;
  725 +}
  726 +
  727 +
  728 +/***********************************************************************//**
  729 + * @brief Save regions into DS9 region file
  730 + *
  731 + * @param[in] filename DS9 region filename.
  732 + *
  733 + * @exception GException::file_open_error
  734 + * File could not be opened.
  735 + *
  736 + * Saves all regions in the container into a DS9 region file.
  737 + ***************************************************************************/
  738 +void GSkyRegions::save(const std::string& filename) const
  739 +{
  740 + // Open file
  741 + std::ofstream ds9file;
  742 + ds9file.open(filename.c_str());
  743 +
  744 + // If file opened correctly, then save regions
  745 + if (ds9file.is_open()) {
  746 +
  747 + // Write global definition
  748 + std::string fileline;
  749 + fileline.append("# Region file format: DS9 version 4.1\n");
  750 + fileline.append("global color=green dashlist=8 3 width=1");
  751 +// fileline.append("font=\"helvetica 10 normal\" select=1");
  752 +// "highlite=1 dash=0 fixed=0 edit=1 move=1"+
  753 +// "delete=1 include=1 source=1\n";
  754 + ds9file << fileline;
  755 +
  756 + // Loop over regions in container
  757 + for (int i = 0; i < size(); ++i) {
  758 + ds9file << m_regions[i]->write() << "\n";
  759 + }
  760 +
  761 + // Close file
  762 + ds9file.close();
  763 +
  764 + // Store filename
  765 + m_filename = filename;
  766 +
  767 + }
  768 +
  769 + // ... otherwise, if file could not be opened then throw an exception
  770 + else {
  771 + throw GException::file_open_error(G_SAVE, filename);
  772 + }
  773 +
  774 + // Return
  775 + return;
  776 +}
  777 +
  778 +
  779 +/***********************************************************************//**
  780 + * @brief Print regions
  781 + *
  782 + * @param[in] chatter Chattiness (defaults to NORMAL).
  783 + * @return String containing region container information.
  784 + *
  785 + * Prints all regions into a string.
  786 + ***************************************************************************/
  787 +std::string GSkyRegions::print(const GChatter& chatter) const
  788 +{
  789 + // Initialise result string
  790 + std::string result;
  791 +
  792 + // Continue only if chatter is not silent
  793 + if (chatter != SILENT) {
  794 +
  795 + // Append header
  796 + result.append("=== GSkyRegions ===");
  797 +
  798 + // Append information
  799 + result.append("\n"+gammalib::parformat("Number of regions"));
  800 + result.append(gammalib::str(size()));
  801 +
  802 + // Append regions
  803 + for (int i = 0; i < size(); ++i) {
  804 + result.append("\n"+m_regions[i]->print(chatter));
  805 + }
  806 +
  807 + } // endif: chatter was not silent
  808 +
  809 + // Return result
  810 + return result;
  811 +}
  812 +
  813 +
  814 +/***********************************************************************//**
  815 + * @brief Tells if direction is contained in one of the regions
  816 + *
  817 + * @param[in] dir A sky direction
  818 + * @return True or False
  819 + *
  820 + * Tells if direction is contained in one of the regions
  821 + ***************************************************************************/
  822 + bool GSkyRegions::contains(const GSkyDir& dir) const
  823 + {
  824 + // Initialise return value
  825 + bool dir_is_in = false;
  826 +
  827 + // Loop over regions
  828 + for (int i = 0; i < size(); ++i) {
  829 + dir_is_in=m_regions[i]->contains(dir);
  830 + if (dir_is_in) break;
  831 + }
  832 +
  833 + // Return result
  834 + return dir_is_in;
  835 + }
  836 +
  837 +
  838 +/***********************************************************************//**
  839 + * @brief Tells if region overlaps one of the regions
  840 + *
  841 + * @param[in] reg A sky region
  842 + * @return True or False
  843 + *
  844 + * Tells if region overlaps one of the regions
  845 + ***************************************************************************/
  846 + bool GSkyRegions::overlaps(const GSkyRegion& reg) const
  847 + {
  848 + // Initialise return value
  849 + bool reg_is_in = false;
  850 +
  851 + // Loop over regions
  852 + for (int i = 0; i < size(); ++i) {
  853 + reg_is_in=m_regions[i]->overlaps(reg);
  854 + if (reg_is_in) break;
  855 + }
  856 +
  857 + // Return result
  858 + return reg_is_in;
  859 + }
  860 +
  861 +
  862 +/*==========================================================================
  863 + = =
  864 + = Private methods =
  865 + = =
  866 + ==========================================================================*/
  867 +
  868 +/***********************************************************************//**
  869 + * @brief Initialise class members
  870 + ***************************************************************************/
  871 +void GSkyRegions::init_members(void)
  872 +{
  873 + // Initialise members
  874 + m_filename.clear();
  875 + m_regions.clear();
  876 +
  877 + // Return
  878 + return;
  879 +}
  880 +
  881 +
  882 +/***********************************************************************//**
  883 + * @brief Copy class members
  884 + *
  885 + * @param[in] regions region container.
  886 + *
  887 + * Makes a copy of all class members. All regions are deep copied, and the
  888 + * linear pointer array for parameter access through the GOptimizerPars
  889 + * base class is set.
  890 + ***************************************************************************/
  891 +void GSkyRegions::copy_members(const GSkyRegions& regions)
  892 +{
  893 + // Copy members
  894 + m_filename = regions.m_filename;
  895 +
  896 + // Copy regions
  897 + m_regions.clear();
  898 + for (int i = 0; i < regions.m_regions.size(); ++i) {
  899 + m_regions.push_back((regions.m_regions[i]->clone()));
  900 + }
  901 +
  902 + // Return
  903 + return;
  904 +}
  905 +
  906 +
  907 +/***********************************************************************//**
  908 + * @brief Delete class members
  909 + *
  910 + * Deallocates all regions. The method loops over the region container and
  911 + * deallocates the memory that has been allocated before.
  912 + ***************************************************************************/
  913 +void GSkyRegions::free_members(void)
  914 +{
  915 + // Free regions
  916 + for (int i = 0; i < m_regions.size(); ++i) {
  917 + delete m_regions[i];
  918 + m_regions[i] = NULL;
  919 + }
  920 +
  921 + // Return
  922 + return;
  923 +}
  924 +
  925 +
  926 +/***********************************************************************//**
  927 + * @brief Return region index by name
  928 + *
  929 + * @param[in] name region name.
  930 + * @return region index (-1 if region name has not been found)
  931 + *
  932 + * Returns region index based on the specified @p name. If no region with the
  933 + * specified @p name is found the method returns -1.
  934 + ***************************************************************************/
  935 +int GSkyRegions::get_index(const std::string& name) const
  936 +{
  937 + // Initialise index
  938 + int index = -1;
  939 +
  940 + // Search region with specified name
  941 + for (int i = 0; i < size(); ++i) {
  942 + if (m_regions[i]->name() == name) {
  943 + index = i;
  944 + break;
  945 + }
  946 + }
  947 +
  948 + // Return index
  949 + return index;
  950 +}
... ...
src/sky/Makefile.am
  View file @7656e43
... ... @@ -9,6 +9,9 @@ INCLUDES = -I$(top_srcdir)/include \
9 9 sources = GSkyDir.cpp \
10 10 GSkyPixel.cpp \
11 11 GSkymap.cpp \
  12 + GSkyRegions.cpp \
  13 + GSkyRegion.cpp \
  14 + GSkyRegionCircle.cpp \
12 15 GWcs.cpp \
13 16 GWcsRegistry.cpp \
14 17 GWcslib.cpp \
... ...
src/support/GTools.cpp
  View file @7656e43
... ... @@ -824,3 +824,25 @@ bool gammalib::file_exists(const std::string&amp; filename)
824 824 // Return result
825 825 return result;
826 826 }
  827 +
  828 +/***********************************************************************//**
  829 + * @brief Checks if a substring is in a string .
  830 + *
  831 + * @param[in] str string you want to search in.
  832 + * @param[in] substring string you are looking for in str.
  833 + *
  834 + * Checks if substring is contained in str
  835 + ***************************************************************************/
  836 +bool gammalib::contains(std::string str, std::string substring)
  837 +{
  838 + // Initialise result
  839 + bool result = false;
  840 +
  841 + // checks if substring is in str
  842 + if (str.find(substring) != std::string::npos){
  843 + result = true;
  844 + }
  845 +
  846 + // Return result
  847 + return result;
  848 +}
... ...
src/xspec/GArf.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GArf.cpp - XSPEC Auxiliary Response File class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GArf.cpp
  23 + * @brief XSPEC Auxiliary Response File class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GArf.hpp"
  32 +#include "GException.hpp"
  33 +#include "GTools.hpp"
  34 +#include "GFitsBinTable.hpp"
  35 +#include "GFitsTableShortCol.hpp"
  36 +#include "GFitsTableFloatCol.hpp"
  37 +
  38 +/* __ Method name definitions ____________________________________________ */
  39 +#define G_AT "GArf::at(int&)"
  40 +
  41 +/* __ Macros _____________________________________________________________ */
  42 +
  43 +/* __ Coding definitions _________________________________________________ */
  44 +
  45 +/* __ Debug definitions __________________________________________________ */
  46 +
  47 +
  48 +/*==========================================================================
  49 + = =
  50 + = Constructors/destructors =
  51 + = =
  52 + ==========================================================================*/
  53 +
  54 +/***********************************************************************//**
  55 + * @brief Void constructor
  56 + ***************************************************************************/
  57 +GArf::GArf(void)
  58 +{
  59 + // Initialise members
  60 + init_members();
  61 +
  62 + // Return
  63 + return;
  64 +}
  65 +
  66 +
  67 +/***********************************************************************//**
  68 + * @brief File constructor
  69 + *
  70 + * @param[in] filename File name.
  71 + ***************************************************************************/
  72 +GArf::GArf(const std::string& filename)
  73 +{
  74 + // Initialise members
  75 + init_members();
  76 +
  77 + // Load ARF file
  78 + load(filename);
  79 +
  80 + // Return
  81 + return;
  82 +}
  83 +
  84 +
  85 +/***********************************************************************//**
  86 + * @brief Energy boundary constructor
  87 + *
  88 + * @param[in] ebds Energy boundaries.
  89 + ***************************************************************************/
  90 +GArf::GArf(const GEbounds& ebds)
  91 +{
  92 + // Initialise members
  93 + init_members();
  94 +
  95 + // Set energy boundaries
  96 + m_ebounds = ebds;
  97 +
  98 + // Initialize spectral response
  99 + m_specresp.assign(ebds.size(), 0.0);
  100 +
  101 + // Return
  102 + return;
  103 +}
  104 +
  105 +
  106 +/***********************************************************************//**
  107 + * @brief Copy constructor
  108 + *
  109 + * @param[in] arf Auxiliary Response File.
  110 + ***************************************************************************/
  111 +GArf::GArf(const GArf& arf)
  112 +{
  113 + // Initialise members
  114 + init_members();
  115 +
  116 + // Copy members
  117 + copy_members(arf);
  118 +
  119 + // Return
  120 + return;
  121 +}
  122 +
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Destructor
  126 + ***************************************************************************/
  127 +GArf::~GArf(void)
  128 +{
  129 + // Free members
  130 + free_members();
  131 +
  132 + // Return
  133 + return;
  134 +}
  135 +
  136 +
  137 +/*==========================================================================
  138 + = =
  139 + = Operators =
  140 + = =
  141 + ==========================================================================*/
  142 +
  143 +/***********************************************************************//**
  144 + * @brief Assignment operator
  145 + *
  146 + * @param[in] arf Auxiliary Response File.
  147 + * @return Auxiliary Response File.
  148 + ***************************************************************************/
  149 +GArf& GArf::operator=(const GArf& arf)
  150 +{
  151 + // Execute only if object is not identical
  152 + if (this != &arf) {
  153 +
  154 + // Free members
  155 + free_members();
  156 +
  157 + // Initialise members
  158 + init_members();
  159 +
  160 + // Copy members
  161 + copy_members(arf);
  162 +
  163 + } // endif: object was not identical
  164 +
  165 + // Return
  166 + return *this;
  167 +}
  168 +
  169 +
  170 +/*==========================================================================
  171 + = =
  172 + = Public methods =
  173 + = =
  174 + ==========================================================================*/
  175 +
  176 +/***********************************************************************//**
  177 + * @brief Clear object.
  178 + *
  179 + * Reset object to a clean initial state.
  180 + ***************************************************************************/
  181 +void GArf::clear(void)
  182 +{
  183 + // Free memory and initialise members
  184 + free_members();
  185 + init_members();
  186 +
  187 + // Return
  188 + return;
  189 +}
  190 +
  191 +
  192 +/***********************************************************************//**
  193 + * @brief Clone object
  194 + *
  195 + * @return Pointer to Auxiliary Response File.
  196 + ***************************************************************************/
  197 +GArf* GArf::clone(void) const
  198 +{
  199 + // Clone client
  200 + return new GArf(*this);
  201 +}
  202 +
  203 +
  204 +/***********************************************************************//**
  205 + * @brief Return content of spectral bin
  206 + *
  207 + * @param[in] index Bin index [0,...,size()-1].
  208 + *
  209 + * @exception GException::out_of_range
  210 + * Spectral bin index is out of range.
  211 + *
  212 + * Returns reference to content of spectral bin with specified @p index.
  213 + ***************************************************************************/
  214 +double& GArf::at(const int& index)
  215 +{
  216 + // Raise exception if index is out of range
  217 + if (index < 0 || index >= size()) {
  218 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  219 + }
  220 +
  221 + // Return reference
  222 + return (m_specresp[index]);
  223 +}
  224 +
  225 +
  226 +/***********************************************************************//**
  227 + * @brief Return content of spectral bin (const version)
  228 + *
  229 + * @param[in] index Bin index [0,...,size()-1].
  230 + *
  231 + * @exception GException::out_of_range
  232 + * Spectral bin index is out of range.
  233 + *
  234 + * Returns reference to content of spectral bin with specified @p index.
  235 + ***************************************************************************/
  236 +const double& GArf::at(const int& index) const
  237 +{
  238 + // Raise exception if index is out of range
  239 + if (index < 0 || index >= size()) {
  240 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  241 + }
  242 +
  243 + // Return reference
  244 + return (m_specresp[index]);
  245 +}
  246 +
  247 +
  248 +/***********************************************************************//**
  249 + * @brief Load Auxiliary Response File
  250 + *
  251 + * @param[in] filename File name.
  252 + ***************************************************************************/
  253 +void GArf::load(const std::string& filename)
  254 +{
  255 + // Clear any existing models
  256 + clear();
  257 +
  258 + // Open FITS file
  259 + GFits file(filename);
  260 +
  261 + // Get ARF table
  262 + GFitsTable* table = file.table("SPECRESP");
  263 +
  264 + // Read ARF data
  265 + read(table);
  266 +
  267 + // Close FITS file
  268 + file.close();
  269 +
  270 + // Store filename
  271 + m_filename = filename;
  272 +
  273 + // Return
  274 + return;
  275 +}
  276 +
  277 +
  278 +/***********************************************************************//**
  279 + * @brief Save Auxiliary Response File
  280 + *
  281 + * @param[in] filename File name.
  282 + * @param[in] clobber Overwrite existing file? (defaults to true)
  283 + ***************************************************************************/
  284 +void GArf::save(const std::string& filename, const bool& clobber) const
  285 +{
  286 + // Open FITS file
  287 + GFits fits;
  288 +
  289 + // Write ARF into file
  290 + write(fits);
  291 +
  292 + // Close FITS file
  293 + fits.saveto(filename, clobber);
  294 +
  295 + // Store filename
  296 + m_filename = filename;
  297 +
  298 + // Return
  299 + return;
  300 +}
  301 +
  302 +
  303 +/***********************************************************************//**
  304 + * @brief Read Auxiliary Response File
  305 + *
  306 + * @param[in] hdu ARF FITS table.
  307 + ***************************************************************************/
  308 +void GArf::read(const GFitsTable* hdu)
  309 +{
  310 + // Continue only if HDU is valid
  311 + if (hdu != NULL) {
  312 +
  313 + // Get pointer to data columns
  314 + const GFitsTableCol* energy_lo = &(*hdu)["ENERG_LO"];
  315 + const GFitsTableCol* energy_hi = &(*hdu)["ENERG_HI"];
  316 + const GFitsTableCol* specresp = &(*hdu)["SPECRESP"];
  317 +
  318 + // Determine effective area conversion factor. Internal
  319 + // units are cm^2
  320 + std::string u_specresp = gammalib::tolower(gammalib::strip_whitespace(specresp->unit()));
  321 + double c_specresp = 1.0;
  322 + if (u_specresp == "m^2" || u_specresp == "m2") {
  323 + c_specresp = 10000.0;
  324 + }
  325 +
  326 + // Extract number of energy bins
  327 + int num = energy_lo->length();
  328 +
  329 + // Set energy bins
  330 + for (int i = 0; i < num; ++i) {
  331 +
  332 + // Append energy bin
  333 + GEnergy emin(energy_lo->real(i), energy_lo->unit());
  334 + GEnergy emax(energy_hi->real(i), energy_hi->unit());
  335 + m_ebounds.append(emin, emax);
  336 +
  337 + // Append effective area value
  338 + double aeff = specresp->real(i) * c_specresp;
  339 + m_specresp.push_back(aeff);
  340 +
  341 + } // endfor: looped over energy bins
  342 +
  343 + } // endif: HDU was valid
  344 +
  345 + // Return
  346 + return;
  347 +}
  348 +
  349 +
  350 +/***********************************************************************//**
  351 + * @brief Write Auxiliary Response File
  352 + *
  353 + * @param[in] fits FITS file.
  354 + ***************************************************************************/
  355 +void GArf::write(GFits& fits) const
  356 +{
  357 + // Set column length
  358 + int length = size();
  359 +
  360 + // Continue only if there are bins
  361 + if (length > 0) {
  362 +
  363 + // Create new binary table
  364 + GFitsBinTable* hdu = new GFitsBinTable;
  365 +
  366 + // Allocate floating point vector columns
  367 + GFitsTableFloatCol energy_lo("ENERG_LO", length);
  368 + GFitsTableFloatCol energy_hi("ENERG_HI", length);
  369 + GFitsTableFloatCol specresp("SPECRESP", length);
  370 +
  371 + // Fill columns
  372 + for (int i = 0; i < length; ++i) {
  373 + energy_lo(i) = m_ebounds.emin(i).keV();
  374 + energy_hi(i) = m_ebounds.emax(i).keV();
  375 + specresp(i) = m_specresp[i];
  376 + }
  377 +
  378 + // Set column units
  379 + energy_lo.unit("keV");
  380 + energy_hi.unit("keV");
  381 + specresp.unit("cm^2");
  382 +
  383 + // Set table attributes
  384 + hdu->extname("SPECRESP");
  385 +
  386 + // Append columns to table
  387 + hdu->append_column(energy_lo);
  388 + hdu->append_column(energy_hi);
  389 + hdu->append_column(specresp);
  390 +
  391 + // Append HDU to FITS file
  392 + fits.append(*hdu);
  393 +
  394 + // Free binary table
  395 + delete hdu;
  396 +
  397 + } // endif: there were data to write
  398 +
  399 + // Return
  400 + return;
  401 +}
  402 +
  403 +
  404 +/***********************************************************************//**
  405 + * @brief Print Auxiliary Response File
  406 + *
  407 + * @param[in] chatter Chattiness (defaults to NORMAL).
  408 + * @return String containing Auxiliary Response File information.
  409 + ***************************************************************************/
  410 +std::string GArf::print(const GChatter& chatter) const
  411 +{
  412 + // Initialise result string
  413 + std::string result;
  414 +
  415 + // Continue only if chatter is not silent
  416 + if (chatter != SILENT) {
  417 +
  418 + // Append header
  419 + result.append("=== GArf ===");
  420 +
  421 + // Append energy boundary information
  422 + result.append("\n"+gammalib::parformat("Number of bins"));
  423 + result.append(gammalib::str(m_ebounds.size()));
  424 + result.append("\n"+gammalib::parformat("Energy range"));
  425 + result.append(m_ebounds.emin().print());
  426 + result.append(" - ");
  427 + result.append(m_ebounds.emax().print());
  428 +
  429 + } // endif: chatter was not silent
  430 +
  431 + // Return result
  432 + return result;
  433 +}
  434 +
  435 +
  436 +/*==========================================================================
  437 + = =
  438 + = Private methods =
  439 + = =
  440 + ==========================================================================*/
  441 +
  442 +/***********************************************************************//**
  443 + * @brief Initialise class members
  444 + ***************************************************************************/
  445 +void GArf::init_members(void)
  446 +{
  447 + // Initialise members
  448 + m_filename.clear();
  449 + m_ebounds.clear();
  450 + m_specresp.clear();
  451 +
  452 + // Return
  453 + return;
  454 +}
  455 +
  456 +
  457 +/***********************************************************************//**
  458 + * @brief Copy class members
  459 + *
  460 + * @param[in] arf Auxiliary Response File.
  461 + ***************************************************************************/
  462 +void GArf::copy_members(const GArf& arf)
  463 +{
  464 + // Copy members
  465 + m_filename = arf.m_filename;
  466 + m_ebounds = arf.m_ebounds;
  467 + m_specresp = arf.m_specresp;
  468 +
  469 + // Return
  470 + return;
  471 +}
  472 +
  473 +
  474 +/***********************************************************************//**
  475 + * @brief Delete class members
  476 + ***************************************************************************/
  477 +void GArf::free_members(void)
  478 +{
  479 + // Return
  480 + return;
  481 +}
... ...
src/xspec/GPha.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GPha.cpp - XSPEC Pulse Height Analyzer class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GPha.cpp
  23 + * @brief XSPEC Pulse Height Analyzer class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GPha.hpp"
  32 +#include "GException.hpp"
  33 +#include "GTools.hpp"
  34 +#include "GFitsBinTable.hpp"
  35 +#include "GFitsTableShortCol.hpp"
  36 +#include "GFitsTableFloatCol.hpp"
  37 +
  38 +/* __ Method name definitions ____________________________________________ */
  39 +#define G_AT "GPha::at(int&)"
  40 +
  41 +/* __ Macros _____________________________________________________________ */
  42 +
  43 +/* __ Coding definitions _________________________________________________ */
  44 +
  45 +/* __ Debug definitions __________________________________________________ */
  46 +
  47 +
  48 +/*==========================================================================
  49 + = =
  50 + = Constructors/destructors =
  51 + = =
  52 + ==========================================================================*/
  53 +
  54 +/***********************************************************************//**
  55 + * @brief Void constructor
  56 + ***************************************************************************/
  57 +GPha::GPha(void)
  58 +{
  59 + // Initialise members
  60 + init_members();
  61 +
  62 + // Return
  63 + return;
  64 +}
  65 +
  66 +
  67 +/***********************************************************************//**
  68 + * @brief File constructor
  69 + *
  70 + * @param[in] filename File name.
  71 + ***************************************************************************/
  72 +GPha::GPha(const std::string& filename)
  73 +{
  74 + // Initialise members
  75 + init_members();
  76 +
  77 + // Load PHA file
  78 + load(filename);
  79 +
  80 + // Return
  81 + return;
  82 +}
  83 +
  84 +
  85 +/***********************************************************************//**
  86 + * @brief Energy boundary constructor
  87 + *
  88 + * @param[in] ebds Energy boundaries.
  89 + ***************************************************************************/
  90 +GPha::GPha(const GEbounds& ebds)
  91 +{
  92 + // Initialise members
  93 + init_members();
  94 +
  95 + // Set energy boundaries
  96 + m_ebounds = ebds;
  97 +
  98 + // Initialize spectrum
  99 + m_counts.assign(ebds.size(), 0.0);
  100 +
  101 + // Return
  102 + return;
  103 +}
  104 +
  105 +
  106 +/***********************************************************************//**
  107 + * @brief Copy constructor
  108 + *
  109 + * @param[in] pha Pulse Height Analyzer spectrum.
  110 + ***************************************************************************/
  111 +GPha::GPha(const GPha& pha)
  112 +{
  113 + // Initialise members
  114 + init_members();
  115 +
  116 + // Copy members
  117 + copy_members(pha);
  118 +
  119 + // Return
  120 + return;
  121 +}
  122 +
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Destructor
  126 + ***************************************************************************/
  127 +GPha::~GPha(void)
  128 +{
  129 + // Free members
  130 + free_members();
  131 +
  132 + // Return
  133 + return;
  134 +}
  135 +
  136 +
  137 +/*==========================================================================
  138 + = =
  139 + = Operators =
  140 + = =
  141 + ==========================================================================*/
  142 +
  143 +/***********************************************************************//**
  144 + * @brief Assignment operator
  145 + *
  146 + * @param[in] pha Pulse Height Analyzer spectrum.
  147 + * @return Pulse Height Analyzer spectrum.
  148 + ***************************************************************************/
  149 +GPha& GPha::operator= (const GPha& pha)
  150 +{
  151 + // Execute only if object is not identical
  152 + if (this != &pha) {
  153 +
  154 + // Free members
  155 + free_members();
  156 +
  157 + // Initialise members
  158 + init_members();
  159 +
  160 + // Copy members
  161 + copy_members(pha);
  162 +
  163 + } // endif: object was not identical
  164 +
  165 + // Return
  166 + return *this;
  167 +}
  168 +
  169 +
  170 +/*==========================================================================
  171 + = =
  172 + = Public methods =
  173 + = =
  174 + ==========================================================================*/
  175 +
  176 +/***********************************************************************//**
  177 + * @brief Clear object.
  178 + *
  179 + * Reset object to a clean initial state.
  180 + ***************************************************************************/
  181 +void GPha::clear(void)
  182 +{
  183 + // Free memory and initialise members
  184 + free_members();
  185 + init_members();
  186 +
  187 + // Return
  188 + return;
  189 +}
  190 +
  191 +
  192 +/***********************************************************************//**
  193 + * @brief Clone object
  194 + * @return Pulse Height Analyzer spectrum.
  195 + ***************************************************************************/
  196 +GPha* GPha::clone(void) const
  197 +{
  198 + // Clone client
  199 + return new GPha(*this);
  200 +}
  201 +
  202 +
  203 +/***********************************************************************//**
  204 + * @brief Return content of spectral bin
  205 + *
  206 + * @param[in] index Bin index [0,...,size()-1].
  207 + *
  208 + * @exception GException::out_of_range
  209 + * Spectral bin index is out of range.
  210 + *
  211 + * Returns reference to content of spectral bin with specified @p index.
  212 + ***************************************************************************/
  213 +double& GPha::at(const int& index)
  214 +{
  215 + // Raise exception if index is out of range
  216 + if (index < 0 || index >= size()) {
  217 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  218 + }
  219 +
  220 + // Return reference
  221 + return (m_counts[index]);
  222 +}
  223 +
  224 +
  225 +/***********************************************************************//**
  226 + * @brief Return content of spectral bin (const version)
  227 + *
  228 + * @param[in] index Bin index [0,...,size()-1].
  229 + *
  230 + * @exception GException::out_of_range
  231 + * Spectral bin index is out of range.
  232 + *
  233 + * Returns reference to content of spectral bin with specified @p index.
  234 + ***************************************************************************/
  235 +const double& GPha::at(const int& index) const
  236 +{
  237 + // Raise exception if index is out of range
  238 + if (index < 0 || index >= size()) {
  239 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  240 + }
  241 +
  242 + // Return reference
  243 + return (m_counts[index]);
  244 +}
  245 +
  246 +
  247 +/***********************************************************************//**
  248 + * @brief Return index for a given energy
  249 + *
  250 + * Returns the number of counts in the spectrum.
  251 + ***************************************************************************/
  252 +/*
  253 +int GPha::index(const GEnergy& energy) const
  254 +{
  255 + // Initialise counts
  256 + double counts = 0.0;
  257 +
  258 + // Compute content
  259 + for (int i = 0; i < m_counts.size(); ++i) {
  260 + counts += m_counts[i];
  261 + }
  262 +
  263 + // Return counts
  264 + return counts;
  265 +}
  266 +*/
  267 +
  268 +
  269 +/***********************************************************************//**
  270 + * @brief Number of counts in spectrum
  271 + *
  272 + * Returns the number of counts in the spectrum.
  273 + ***************************************************************************/
  274 +double GPha::counts(void) const
  275 +{
  276 + // Initialise counts
  277 + double counts = 0.0;
  278 +
  279 + // Compute content
  280 + for (int i = 0; i < m_counts.size(); ++i) {
  281 + counts += m_counts[i];
  282 + }
  283 +
  284 + // Return counts
  285 + return counts;
  286 +}
  287 +
  288 +
  289 +/***********************************************************************//**
  290 + * @brief Fill spectrum with a value.
  291 + *
  292 + * @param[in] energy Energy.
  293 + * @param[in] value Fill value (defaults to 1.0).
  294 + *
  295 + * Fills the specified @p value at a given @p energy in the spectrum.
  296 + ***************************************************************************/
  297 +void GPha::fill(const GEnergy& energy, const double& value)
  298 +{
  299 + // Get index
  300 + int index = m_ebounds.index(energy);
  301 +
  302 + // If index is not valid, then check for underflow, overflow or outflow
  303 + if (index == -1) {
  304 + if (energy < m_ebounds.emin()) {
  305 + m_underflow += value;
  306 + }
  307 + else if (energy >= m_ebounds.emax()) {
  308 + m_overflow += value;
  309 + }
  310 + else {
  311 + m_outflow += value;
  312 + }
  313 + }
  314 +
  315 + // ... otherwise fill the histogram
  316 + else {
  317 + m_counts[index] += value;
  318 + }
  319 +
  320 + // Return
  321 + return;
  322 +}
  323 +
  324 +
  325 +/***********************************************************************//**
  326 + * @brief Load Pulse Height Analyzer spectrum
  327 + *
  328 + * @param[in] filename File name.
  329 + ***************************************************************************/
  330 +void GPha::load(const std::string& filename)
  331 +{
  332 + // Clear any existing models
  333 + clear();
  334 +
  335 + // Open FITS file
  336 + GFits file(filename);
  337 +
  338 + // Read energy boundaries
  339 + m_ebounds.load(filename);
  340 +
  341 + // Get PHA table
  342 + GFitsTable* table = file.table("SPECTRUM");
  343 +
  344 + // Read PHA data
  345 + read(table);
  346 +
  347 + // Close FITS file
  348 + file.close();
  349 +
  350 + // Store filename
  351 + m_filename = filename;
  352 +
  353 + // Return
  354 + return;
  355 +}
  356 +
  357 +
  358 +/***********************************************************************//**
  359 + * @brief Save Pulse Height Analyzer spectrum
  360 + *
  361 + * @param[in] filename File name.
  362 + * @param[in] clobber Overwrite existing file? (defaults to true)
  363 + ***************************************************************************/
  364 +void GPha::save(const std::string& filename, const bool& clobber) const
  365 +{
  366 + // Open FITS file
  367 + GFits fits;
  368 +
  369 + // Write PHA into file
  370 + write(fits);
  371 +
  372 + // Close FITS file
  373 + fits.saveto(filename, clobber);
  374 +
  375 + // Store filename
  376 + m_filename = filename;
  377 +
  378 + // Return
  379 + return;
  380 +}
  381 +
  382 +
  383 +/***********************************************************************//**
  384 + * @brief Read Pulse Height Analyzer spectrum
  385 + *
  386 + * @param[in] hdu PHA FITS table.
  387 + ***************************************************************************/
  388 +void GPha::read(const GFitsTable* hdu)
  389 +{
  390 + // Continue only if HDU is valid
  391 + if (hdu != NULL) {
  392 +
  393 + // Set expected column length
  394 + int length = m_ebounds.size();
  395 +
  396 + // Initialize spectrum
  397 + m_counts.assign(length, 0.0);
  398 +
  399 + // Get pointer to data column
  400 + const GFitsTableCol* col_data = &(*hdu)["COUNTS"];
  401 +
  402 + // Check whether column length is okay
  403 + //TODO
  404 +
  405 + // Copy data
  406 + for (int i = 0; i < length; ++i) {
  407 + m_counts[i] = col_data->real(i);
  408 + }
  409 +
  410 + } // endif: HDU was valid
  411 +
  412 + // Return
  413 + return;
  414 +}
  415 +
  416 +
  417 +/***********************************************************************//**
  418 + * @brief Write Pulse Height Analyzer spectrum
  419 + *
  420 + * @param[in] fits FITS file.
  421 + ***************************************************************************/
  422 +void GPha::write(GFits& fits) const
  423 +{
  424 + // Set column length
  425 + int length = size();
  426 +
  427 + // Continue only if there are bins
  428 + if (length > 0) {
  429 +
  430 + // Create new binary table
  431 + GFitsBinTable* hdu = new GFitsBinTable;
  432 +
  433 + // Allocate floating point vector columns
  434 + GFitsTableShortCol col_chan("CHANNEL", length);
  435 + GFitsTableFloatCol col_data("COUNTS", length);
  436 + GFitsTableFloatCol col_stat("STAT_ERR", length);
  437 + GFitsTableFloatCol col_syst("SYS_ERR", length);
  438 + GFitsTableShortCol col_qual("QUALITY", length);
  439 + GFitsTableShortCol col_grpg("GROUPING", length);
  440 + GFitsTableFloatCol col_area("AREASCAL", length);
  441 + GFitsTableFloatCol col_back("BACKSCAL", length);
  442 +
  443 + // Fill columns
  444 + for (int i = 0; i < length; ++i) {
  445 + col_chan(i) = i+1; // Channels start at 1
  446 + col_data(i) = float(m_counts[i]);
  447 + }
  448 +
  449 + // Set table attributes
  450 + hdu->extname("SPECTRUM");
  451 +
  452 + // Append columns to table
  453 + hdu->append_column(col_chan);
  454 + hdu->append_column(col_data);
  455 + hdu->append_column(col_stat);
  456 + hdu->append_column(col_syst);
  457 + hdu->append_column(col_qual);
  458 + hdu->append_column(col_grpg);
  459 + hdu->append_column(col_area);
  460 + hdu->append_column(col_back);
  461 +
  462 + // Append HDU to FITS file
  463 + fits.append(*hdu);
  464 +
  465 + // Free binary table
  466 + delete hdu;
  467 +
  468 + // Append energy boundaries
  469 + m_ebounds.write(&fits);
  470 +
  471 + } // endif: there were data to write
  472 +
  473 + // Return
  474 + return;
  475 +}
  476 +
  477 +
  478 +/***********************************************************************//**
  479 + * @brief Print Pulse Height Analyzer spectrum
  480 + *
  481 + * @param[in] chatter Chattiness (defaults to NORMAL).
  482 + * @return String containing Pulse Height Analyzer spectrum information.
  483 + ***************************************************************************/
  484 +std::string GPha::print(const GChatter& chatter) const
  485 +{
  486 + // Initialise result string
  487 + std::string result;
  488 +
  489 + // Continue only if chatter is not silent
  490 + if (chatter != SILENT) {
  491 +
  492 + // Append header
  493 + result.append("=== GPha ===");
  494 +
  495 + // Append energy boundary information
  496 + result.append("\n"+gammalib::parformat("Number of bins"));
  497 + result.append(gammalib::str(m_ebounds.size()));
  498 + result.append("\n"+gammalib::parformat("Energy range"));
  499 + result.append(m_ebounds.emin().print());
  500 + result.append(" - ");
  501 + result.append(m_ebounds.emax().print());
  502 +
  503 + // Append content information
  504 + result.append("\n"+gammalib::parformat("Total number of counts"));
  505 + result.append(gammalib::str(counts()));
  506 + result.append("\n"+gammalib::parformat("Underflow counts"));
  507 + result.append(gammalib::str(m_underflow));
  508 + result.append("\n"+gammalib::parformat("Overflow counts"));
  509 + result.append(gammalib::str(m_overflow));
  510 + result.append("\n"+gammalib::parformat("Outflow counts"));
  511 + result.append(gammalib::str(m_outflow));
  512 +
  513 + } // endif: chatter was not silent
  514 +
  515 + // Return result
  516 + return result;
  517 +}
  518 +
  519 +
  520 +/*==========================================================================
  521 + = =
  522 + = Private methods =
  523 + = =
  524 + ==========================================================================*/
  525 +
  526 +/***********************************************************************//**
  527 + * @brief Initialise class members
  528 + ***************************************************************************/
  529 +void GPha::init_members(void)
  530 +{
  531 + // Initialise members
  532 + m_filename.clear();
  533 + m_ebounds.clear();
  534 + m_counts.clear();
  535 + m_underflow = 0.0;
  536 + m_overflow = 0.0;
  537 + m_outflow = 0.0;
  538 +
  539 + // Return
  540 + return;
  541 +}
  542 +
  543 +
  544 +/***********************************************************************//**
  545 + * @brief Copy class members
  546 + *
  547 + * @param[in] pha Pulse Height Analyzer spectrum.
  548 + ***************************************************************************/
  549 +void GPha::copy_members(const GPha& pha)
  550 +{
  551 + // Copy members
  552 + m_filename = pha.m_filename;
  553 + m_ebounds = pha.m_ebounds;
  554 + m_counts = pha.m_counts;
  555 + m_underflow = pha.m_underflow;
  556 + m_overflow = pha.m_overflow;
  557 + m_outflow = pha.m_outflow;
  558 +
  559 + // Return
  560 + return;
  561 +}
  562 +
  563 +
  564 +/***********************************************************************//**
  565 + * @brief Delete class members
  566 + ***************************************************************************/
  567 +void GPha::free_members(void)
  568 +{
  569 + // Return
  570 + return;
  571 +}
... ...
src/xspec/GRmf.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * GRmf.cpp - XSPEC Redistribution Matrix File class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GRmf.cpp
  23 + * @brief XSPEC Redistribution Matrix File class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GRmf.hpp"
  32 +#include "GException.hpp"
  33 +#include "GTools.hpp"
  34 +#include "GFitsBinTable.hpp"
  35 +#include "GFitsTableShortCol.hpp"
  36 +#include "GFitsTableFloatCol.hpp"
  37 +
  38 +/* __ Method name definitions ____________________________________________ */
  39 +#define G_AT "GRmf::at(int&, int&)"
  40 +
  41 +/* __ Macros _____________________________________________________________ */
  42 +
  43 +/* __ Coding definitions _________________________________________________ */
  44 +
  45 +/* __ Debug definitions __________________________________________________ */
  46 +
  47 +
  48 +/*==========================================================================
  49 + = =
  50 + = Constructors/destructors =
  51 + = =
  52 + ==========================================================================*/
  53 +
  54 +/***********************************************************************//**
  55 + * @brief Void constructor
  56 + ***************************************************************************/
  57 +GRmf::GRmf(void)
  58 +{
  59 + // Initialise members
  60 + init_members();
  61 +
  62 + // Return
  63 + return;
  64 +}
  65 +
  66 +
  67 +/***********************************************************************//**
  68 + * @brief File constructor
  69 + *
  70 + * @param[in] filename File name.
  71 + ***************************************************************************/
  72 +GRmf::GRmf(const std::string& filename)
  73 +{
  74 + // Initialise members
  75 + init_members();
  76 +
  77 + // Load RMF file
  78 + load(filename);
  79 +
  80 + // Return
  81 + return;
  82 +}
  83 +
  84 +
  85 +/***********************************************************************//**
  86 + * @brief Energy boundary constructor
  87 + *
  88 + * @param[in] etrue True energy boundaries.
  89 + * @param[in] emeasured Measured energy boundaries.
  90 + ***************************************************************************/
  91 +GRmf::GRmf(const GEbounds& etrue, const GEbounds& emeasured)
  92 +{
  93 + // Initialise members
  94 + init_members();
  95 +
  96 + // Set energy boundaries
  97 + m_ebds_true = etrue;
  98 + m_ebds_measured = emeasured;
  99 +
  100 + // Initialize matrix
  101 + m_matrix = GMatrixSparse(m_ebds_true.size(), m_ebds_measured.size());
  102 +
  103 + // Return
  104 + return;
  105 +}
  106 +
  107 +
  108 +/***********************************************************************//**
  109 + * @brief Copy constructor
  110 + *
  111 + * @param[in] rmf Redistribution Matrix File.
  112 + ***************************************************************************/
  113 +GRmf::GRmf(const GRmf& rmf)
  114 +{
  115 + // Initialise members
  116 + init_members();
  117 +
  118 + // Copy members
  119 + copy_members(rmf);
  120 +
  121 + // Return
  122 + return;
  123 +}
  124 +
  125 +
  126 +/***********************************************************************//**
  127 + * @brief Destructor
  128 + ***************************************************************************/
  129 +GRmf::~GRmf(void)
  130 +{
  131 + // Free members
  132 + free_members();
  133 +
  134 + // Return
  135 + return;
  136 +}
  137 +
  138 +
  139 +/*==========================================================================
  140 + = =
  141 + = Operators =
  142 + = =
  143 + ==========================================================================*/
  144 +
  145 +/***********************************************************************//**
  146 + * @brief Assignment operator
  147 + *
  148 + * @param[in] rmf Redistribution Matrix File.
  149 + * @return Redistribution Matrix File.
  150 + ***************************************************************************/
  151 +GRmf& GRmf::operator=(const GRmf& rmf)
  152 +{
  153 + // Execute only if object is not identical
  154 + if (this != &rmf) {
  155 +
  156 + // Free members
  157 + free_members();
  158 +
  159 + // Initialise members
  160 + init_members();
  161 +
  162 + // Copy members
  163 + copy_members(rmf);
  164 +
  165 + } // endif: object was not identical
  166 +
  167 + // Return
  168 + return *this;
  169 +}
  170 +
  171 +
  172 +/*==========================================================================
  173 + = =
  174 + = Public methods =
  175 + = =
  176 + ==========================================================================*/
  177 +
  178 +/***********************************************************************//**
  179 + * @brief Clear object.
  180 + *
  181 + * Reset object to a clean initial state.
  182 + ***************************************************************************/
  183 +void GRmf::clear(void)
  184 +{
  185 + // Free memory and initialise members
  186 + free_members();
  187 + init_members();
  188 +
  189 + // Return
  190 + return;
  191 +}
  192 +
  193 +
  194 +/***********************************************************************//**
  195 + * @brief Clone object
  196 + *
  197 + * @return Redistribution Matrix File.
  198 + ***************************************************************************/
  199 +GRmf* GRmf::clone(void) const
  200 +{
  201 + // Clone object
  202 + return new GRmf(*this);
  203 +}
  204 +
  205 +
  206 +/***********************************************************************//**
  207 + * @brief Return content of redistribution matrix bin
  208 + *
  209 + * @param[in] itrue True energy index [0,...,ntrue()-1].
  210 + * @param[in] imeasured Measured energy index [0,...,nmeasured()-1].
  211 + *
  212 + * @exception GException::out_of_range
  213 + * Bin index is out of range.
  214 + *
  215 + * Returns reference to content of redistribution matrix bin bin with true
  216 + * energy index @p itrue and measured energy index @p imeasured.
  217 + ***************************************************************************/
  218 +double& GRmf::at(const int& itrue, const int& imeasured)
  219 +{
  220 + // Raise exception if indices are out of range
  221 + if (itrue < 0 || itrue >= ntrue()) {
  222 + throw GException::out_of_range(G_AT, itrue, 0, ntrue()-1);
  223 + }
  224 + // Raise exception if indices are out of range
  225 + if (imeasured < 0 || imeasured >= nmeasured()) {
  226 + throw GException::out_of_range(G_AT, imeasured, 0, nmeasured()-1);
  227 + }
  228 +
  229 + // Return reference
  230 + return (m_matrix(itrue, imeasured));
  231 +}
  232 +
  233 +
  234 +/***********************************************************************//**
  235 + * @brief Return content of redistribution matrix bin (const version)
  236 + *
  237 + * @param[in] itrue True energy index [0,...,ntrue()-1].
  238 + * @param[in] imeasured Measured energy index [0,...,nmeasured()-1].
  239 + *
  240 + * @exception GException::out_of_range
  241 + * Bin index is out of range.
  242 + *
  243 + * Returns reference to content of redistribution matrix bin bin with true
  244 + * energy index @p itrue and measured energy index @p imeasured.
  245 + ***************************************************************************/
  246 +const double& GRmf::at(const int& itrue, const int& imeasured) const
  247 +{
  248 + // Raise exception if indices are out of range
  249 + if (itrue < 0 || itrue >= ntrue()) {
  250 + throw GException::out_of_range(G_AT, itrue, 0, ntrue()-1);
  251 + }
  252 + // Raise exception if indices are out of range
  253 + if (imeasured < 0 || imeasured >= nmeasured()) {
  254 + throw GException::out_of_range(G_AT, imeasured, 0, nmeasured()-1);
  255 + }
  256 +
  257 + // Return reference
  258 + return (m_matrix(itrue, imeasured));
  259 +}
  260 +
  261 +
  262 +/***********************************************************************//**
  263 + * @brief Load Redistribution Matrix File
  264 + *
  265 + * @param[in] filename File name.
  266 + ***************************************************************************/
  267 +void GRmf::load(const std::string& filename)
  268 +{
  269 + // Clear any existing models
  270 + clear();
  271 +
  272 + // Open FITS file
  273 + GFits file(filename);
  274 +
  275 + // Read measured energy boundaries
  276 + m_ebds_measured.load(filename);
  277 +
  278 + // Get RMF table
  279 + GFitsTable* table = file.table("MATRIX");
  280 +
  281 + // Read RMF data
  282 + read(table);
  283 +
  284 + // Close FITS file
  285 + file.close();
  286 +
  287 + // Store filename
  288 + m_filename = filename;
  289 +
  290 + // Return
  291 + return;
  292 +}
  293 +
  294 +
  295 +/***********************************************************************//**
  296 + * @brief Save Redistribution Matrix File
  297 + *
  298 + * @param[in] filename File name.
  299 + * @param[in] clobber Overwrite existing file? (defaults to true)
  300 + ***************************************************************************/
  301 +void GRmf::save(const std::string& filename, const bool& clobber) const
  302 +{
  303 + // Open FITS file
  304 + GFits fits;
  305 +
  306 + // Write RMF into file
  307 + write(fits);
  308 +
  309 + // Close FITS file
  310 + fits.saveto(filename, clobber);
  311 +
  312 + // Store filename
  313 + m_filename = filename;
  314 +
  315 + // Return
  316 + return;
  317 +}
  318 +
  319 +
  320 +/***********************************************************************//**
  321 + * @brief Read Redistribution Matrix File
  322 + *
  323 + * @param[in] hdu RMF FITS table.
  324 + ***************************************************************************/
  325 +void GRmf::read(const GFitsTable* hdu)
  326 +{
  327 + // Continue only if HDU is valid
  328 + if (hdu != NULL) {
  329 +
  330 + // Initialise members
  331 + m_ebds_true.clear();
  332 +
  333 + // Get pointer to data columns
  334 + const GFitsTableCol* energy_lo = &(*hdu)["ENERG_LO"];
  335 + const GFitsTableCol* energy_hi = &(*hdu)["ENERG_HI"];
  336 + const GFitsTableCol* n_grp = &(*hdu)["N_GRP"];
  337 + const GFitsTableCol* f_chan = &(*hdu)["F_CHAN"];
  338 + const GFitsTableCol* n_chan = &(*hdu)["N_CHAN"];
  339 + const GFitsTableCol* matrix = &(*hdu)["MATRIX"];
  340 +
  341 + // Set matrix rows and columns
  342 + int rows = energy_lo->length();
  343 + int columns = m_ebds_measured.size();
  344 +
  345 + // Initialize matrix
  346 + m_matrix = GMatrixSparse(rows, columns);
  347 +
  348 + // Set true energy bins
  349 + for (int itrue = 0; itrue < rows; ++itrue) {
  350 +
  351 + // Append energy bin
  352 + GEnergy emin(energy_lo->real(itrue), energy_lo->unit());
  353 + GEnergy emax(energy_hi->real(itrue), energy_hi->unit());
  354 + m_ebds_true.append(emin, emax);
  355 +
  356 + // Loop over groups
  357 + int icolumn = 0;
  358 + int ngroups = n_grp->integer(itrue);
  359 + for (int igroup = 0; igroup < ngroups; ++igroup) {
  360 +
  361 + // Get start column index and number of columns
  362 + int imeasured = f_chan->integer(itrue, igroup);
  363 + int nvalues = n_chan->integer(itrue, igroup);
  364 +
  365 + // Get values
  366 + for (int i = 0; i < nvalues; ++i, ++imeasured, ++icolumn) {
  367 + m_matrix(itrue, imeasured) = matrix->real(itrue, icolumn);
  368 + }
  369 +
  370 + } // endfor: looped over groups
  371 +
  372 + } // endfor: looped over true energy bins
  373 +
  374 + } // endif: HDU was valid
  375 +
  376 + // Return
  377 + return;
  378 +}
  379 +
  380 +
  381 +/***********************************************************************//**
  382 + * @brief Write Redistribution Matrix File
  383 + *
  384 + * @param[in] fits FITS file.
  385 + ***************************************************************************/
  386 +void GRmf::write(GFits& fits) const
  387 +{
  388 + // Set table length
  389 + int length = ntrue();
  390 +
  391 + // Continue only if there are true energies
  392 + if (length > 0) {
  393 +
  394 + // Compute maximum number of matrix elements
  395 + int maxcolumns = 1;
  396 + for (int i = 0; i < length; i++) {
  397 +
  398 + // Search first non-zero element
  399 + int ifirst = 0;
  400 + for (; ifirst < nmeasured(); ++ifirst) {
  401 + if (m_matrix(i, ifirst) != 0.0) {
  402 + break;
  403 + }
  404 + }
  405 +
  406 + // Search last non-zero element
  407 + int ilast = nmeasured() - 1;
  408 + for (; ilast >= 0; --ilast) {
  409 + if (m_matrix(i, ilast) != 0.0) {
  410 + break;
  411 + }
  412 + }
  413 +
  414 + // Compute number of non-zero matrix elements
  415 + int num = ilast - ifirst + 1;
  416 +
  417 + // Set maximum
  418 + if (num > maxcolumns) {
  419 + maxcolumns = num;
  420 + }
  421 +
  422 + } // endfor: looped over true energies
  423 +
  424 + // Create new binary table
  425 + GFitsBinTable* hdu = new GFitsBinTable;
  426 +
  427 + // Allocate columns
  428 + GFitsTableFloatCol energy_lo("ENERG_LO", length);
  429 + GFitsTableFloatCol energy_hi("ENERG_HI", length);
  430 + GFitsTableShortCol n_grp("N_GRP", length);
  431 + GFitsTableShortCol f_chan("F_CHAN", length);
  432 + GFitsTableShortCol n_chan("N_CHAN", length);
  433 + GFitsTableFloatCol matrix("MATRIX", length, maxcolumns);
  434 +
  435 + // Fill columns
  436 + for (int itrue = 0; itrue < length; ++itrue) {
  437 +
  438 + // Search first non-zero element
  439 + int ifirst = 0;
  440 + for (; ifirst < nmeasured(); ++ifirst) {
  441 + if (m_matrix(itrue, ifirst) != 0.0) {
  442 + break;
  443 + }
  444 + }
  445 +
  446 + // Search last non-zero element
  447 + int ilast = nmeasured() - 1;
  448 + for (; ilast >= 0; --ilast) {
  449 + if (m_matrix(itrue, ilast) != 0.0) {
  450 + break;
  451 + }
  452 + }
  453 +
  454 + // Compute number of non-zero matrix elements
  455 + int num = ilast - ifirst + 1;
  456 + if (num < 0) {
  457 + ifirst = 0;
  458 + num = 0;
  459 + }
  460 +
  461 + // Set energy bin
  462 + energy_lo(itrue) = m_ebds_true.emin(itrue).keV();
  463 + energy_hi(itrue) = m_ebds_true.emax(itrue).keV();
  464 +
  465 + // Set group information
  466 + n_grp(itrue) = 1;
  467 + f_chan(itrue) = ifirst;
  468 + n_chan(itrue) = num;
  469 +
  470 + // Set matrix elements
  471 + for (int i = 0; i < num; ++i, ++ifirst) {
  472 + matrix(itrue, i) = m_matrix(itrue, ifirst);
  473 + }
  474 +
  475 + } // endfor: looped over true energy bins
  476 +
  477 + // Set column units
  478 + energy_lo.unit("keV");
  479 + energy_hi.unit("keV");
  480 +
  481 + // Set table attributes
  482 + hdu->extname("MATRIX");
  483 +
  484 + // Append columns to table
  485 + hdu->append_column(energy_lo);
  486 + hdu->append_column(energy_hi);
  487 + hdu->append_column(n_grp);
  488 + hdu->append_column(f_chan);
  489 + hdu->append_column(n_chan);
  490 + hdu->append_column(matrix);
  491 +
  492 + // Append HDU to FITS file
  493 + fits.append(*hdu);
  494 +
  495 + // Free binary table
  496 + delete hdu;
  497 +
  498 + // Append measured energy boundaries
  499 + m_ebds_measured.write(&fits);
  500 +
  501 + } // endif: there were data to write
  502 +
  503 + // Return
  504 + return;
  505 +}
  506 +
  507 +
  508 +/***********************************************************************//**
  509 + * @brief Print Redistribution Matrix File
  510 + *
  511 + * @param[in] chatter Chattiness (defaults to NORMAL).
  512 + * @return String containing Redistribution Matrix File information.
  513 + ***************************************************************************/
  514 +std::string GRmf::print(const GChatter& chatter) const
  515 +{
  516 + // Initialise result string
  517 + std::string result;
  518 +
  519 + // Continue only if chatter is not silent
  520 + if (chatter != SILENT) {
  521 +
  522 + // Append header
  523 + result.append("=== GRmf ===");
  524 +
  525 + // Append energy boundary information
  526 + result.append("\n"+gammalib::parformat("Number of true energy bins"));
  527 + result.append(gammalib::str(m_ebds_true.size()));
  528 + result.append("\n"+gammalib::parformat("Number of measured bins"));
  529 + result.append(gammalib::str(m_ebds_measured.size()));
  530 + result.append("\n"+gammalib::parformat("True energy range"));
  531 + result.append(m_ebds_true.emin().print());
  532 + result.append(" - ");
  533 + result.append(m_ebds_true.emax().print());
  534 + result.append("\n"+gammalib::parformat("Measured energy range"));
  535 + result.append(m_ebds_measured.emin().print());
  536 + result.append(" - ");
  537 + result.append(m_ebds_measured.emax().print());
  538 +
  539 + } // endif: chatter was not silent
  540 +
  541 + // Return result
  542 + return result;
  543 +}
  544 +
  545 +
  546 +/*==========================================================================
  547 + = =
  548 + = Private methods =
  549 + = =
  550 + ==========================================================================*/
  551 +
  552 +/***********************************************************************//**
  553 + * @brief Initialise class members
  554 + ***************************************************************************/
  555 +void GRmf::init_members(void)
  556 +{
  557 + // Initialise members
  558 + m_filename.clear();
  559 + m_ebds_true.clear();
  560 + m_ebds_measured.clear();
  561 + m_matrix.clear();
  562 +
  563 + // Return
  564 + return;
  565 +}
  566 +
  567 +
  568 +/***********************************************************************//**
  569 + * @brief Copy class members
  570 + *
  571 + * @param[in] rmf Redistribution Matrix File.
  572 + ***************************************************************************/
  573 +void GRmf::copy_members(const GRmf& rmf)
  574 +{
  575 + // Copy members
  576 + m_filename = rmf.m_filename;
  577 + m_ebds_true = rmf.m_ebds_true;
  578 + m_ebds_measured = rmf.m_ebds_measured;
  579 + m_matrix = rmf.m_matrix;
  580 +
  581 + // Return
  582 + return;
  583 +}
  584 +
  585 +
  586 +/***********************************************************************//**
  587 + * @brief Delete class members
  588 + ***************************************************************************/
  589 +void GRmf::free_members(void)
  590 +{
  591 + // Return
  592 + return;
  593 +}
... ...
src/xspec/Makefile.am 0 → 100644
  View file @7656e43
  1 +## Process this file with automake to produce Makefile.in
  2 +
  3 +# Set include directories
  4 +INCLUDES = -I$(top_srcdir)/include \
  5 + -I$(top_builddir)/src \
  6 + -I$(top_srcdir)/src/support
  7 +
  8 +# Define sources for this directory
  9 +sources = GArf.cpp \
  10 + GPha.cpp \
  11 + GRmf.cpp
  12 +
  13 +# Build libtool library
  14 +noinst_LTLIBRARIES = libxspec.la
  15 +libxspec_la_SOURCES = $(sources)
  16 +libxspec_la_LDFLAGS = -version-info $(GAMMALIB_LT_VERSION)
... ...
test/Makefile.am
  View file @7656e43
... ... @@ -71,6 +71,7 @@ check_PROGRAMS = test_GSupport \
71 71 test_GFits \
72 72 test_GXml \
73 73 test_GVO \
  74 + test_GXspec \
74 75 test_GApplication \
75 76 test_GModel \
76 77 test_GSky \
... ... @@ -92,6 +93,7 @@ TESTS = test_GSupport \
92 93 test_GFits \
93 94 test_GXml \
94 95 test_GVO \
  96 + test_GXspec \
95 97 test_GApplication \
96 98 test_GModel \
97 99 test_GSky \
... ... @@ -146,6 +148,11 @@ test_GVO_LDFLAGS = @LDFLAGS@
146 148 test_GVO_CPPFLAGS = @CPPFLAGS@
147 149 test_GVO_LDADD = $(top_srcdir)/src/libgamma.la
148 150 #
  151 +test_GXspec_SOURCES = test_GXspec.cpp
  152 +test_GXspec_LDFLAGS = @LDFLAGS@
  153 +test_GXspec_CPPFLAGS = @CPPFLAGS@
  154 +test_GXspec_LDADD = $(top_srcdir)/src/libgamma.la
  155 +#
149 156 test_GApplication_SOURCES = test_GApplication.cpp
150 157 test_GApplication_LDFLAGS = @LDFLAGS@
151 158 test_GApplication_CPPFLAGS = @CPPFLAGS@
... ...
test/data/test_circle_radeczerozero.reg 0 → 100644
  View file @7656e43
  1 +# Region file format: DS9 version 4.1
  2 +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
  3 +fk5
  4 +circle(0.,0.,10.)
... ...
test/test_GSky.cpp
  View file @7656e43
1 1 /***************************************************************************
2 2 * test_GSky.cpp - test GSky classes *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2013 by Jurgen Knodlseder *
  4 + * copyright (C) 2010-2013 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -21,7 +21,7 @@
21 21 /**
22 22 * @file test_GSky.cpp
23 23 * @brief Testing of sky module.
24   - * @author J. Knodlseder
  24 + * @author Juergen Knoedlseder
25 25 */
26 26  
27 27 /* __ Includes ___________________________________________________________ */
... ... @@ -41,9 +41,11 @@
41 41 /* __ Debug definitions __________________________________________________ */
42 42 #define G_WCS_FORTH_BACK_PIXEL_DEBUG
43 43 #define G_WCS_COPY_DEBUG
  44 +
  45 +
44 46 /***********************************************************************//**
45   -* @brief Set parameters and tests
46   -***************************************************************************/
  47 + * @brief Set parameters and tests
  48 + ***************************************************************************/
47 49 void TestGSky::set(void){
48 50  
49 51 // Test name
... ... @@ -55,7 +57,9 @@ void TestGSky::set(void){
55 57 add_test(static_cast<pfunction>(&TestGSky::test_GSkymap_healpix_io),"Test Healpix GSkymap I/O");
56 58 add_test(static_cast<pfunction>(&TestGSky::test_GSkymap_wcs_construct),"Test WCS GSkymap constructors");
57 59 add_test(static_cast<pfunction>(&TestGSky::test_GSkymap_wcs_io),"Test WCS GSkymap I/O");
58   -
  60 + add_test(static_cast<pfunction>(&TestGSky::test_GSkyRegions_io),"Test the GSkyregion class - io");
  61 + add_test(static_cast<pfunction>(&TestGSky::test_GSkyRegionCircle_construct),"Test the GSkyregionCircle class - constuct");
  62 + add_test(static_cast<pfunction>(&TestGSky::test_GSkyRegionCircle_logic),"Test the GSkyregionCircle class - logic ");
59 63 return;
60 64 }
61 65  
... ... @@ -353,7 +357,7 @@ void TestGSky::test_GWcslib(void)
353 357  
354 358  
355 359 /***************************************************************************
356   - * Test: GSkymap_healpix_construct *
  360 + * @brief GSkymap_healpix_construct
357 361 ***************************************************************************/
358 362 void TestGSky::test_GSkymap_healpix_construct(void)
359 363 {
... ... @@ -499,7 +503,7 @@ void TestGSky::test_GSkymap_healpix_construct(void)
499 503  
500 504  
501 505 /***************************************************************************
502   - * Test: GSkymap_healpix_io *
  506 + * @brief GSkymap_healpix_io
503 507 ***************************************************************************/
504 508 void TestGSky::test_GSkymap_healpix_io(void)
505 509 {
... ... @@ -567,7 +571,7 @@ void TestGSky::test_GSkymap_healpix_io(void)
567 571  
568 572  
569 573 /***************************************************************************
570   - * Test: GSkymap_wcs_construct *
  574 + * @brief GSkymap_wcs_construct
571 575 ***************************************************************************/
572 576 void TestGSky::test_GSkymap_wcs_construct(void)
573 577 {
... ... @@ -626,7 +630,7 @@ void TestGSky::test_GSkymap_wcs_construct(void)
626 630  
627 631  
628 632 /***************************************************************************
629   - * Test: GSkymap_wcs_io *
  633 + * @brief GSkymap_wcs_io
630 634 ***************************************************************************/
631 635 void TestGSky::test_GSkymap_wcs_io(void)
632 636 {
... ... @@ -700,7 +704,190 @@ void TestGSky::test_GSkymap_wcs_io(void)
700 704  
701 705  
702 706 /***************************************************************************
703   - * Main test function *
  707 + * @brief GSkyRegionCircle_construct
  708 + ***************************************************************************/
  709 +void TestGSky::test_GSkyRegionCircle_construct(void)
  710 +{
  711 + // Define region for comparison
  712 + GSkyDir refdir_radeczerozero = GSkyDir();
  713 + refdir_radeczerozero.radec_deg(0,0);
  714 + double refradius = 10.;
  715 +
  716 + // Test constructing:
  717 + test_try("Test constructor");
  718 + try {
  719 + GSkyRegionCircle circle(refdir_radeczerozero,refradius);
  720 + test_try_success();
  721 + }
  722 + catch (std::exception &e) {
  723 + test_try_failure(e);
  724 + }
  725 +
  726 + //Test constructing with radius 0
  727 + test_try("Test constructor2");
  728 + try {
  729 + GSkyRegionCircle circle(refdir_radeczerozero,-1);
  730 + test_try_failure();
  731 + }
  732 + catch (GException::invalid_argument &e) {
  733 + test_try_success();
  734 + }
  735 + catch (std::exception &e) {
  736 + test_try_failure(e);
  737 + }
  738 +
  739 + //Test constructing with radius 0
  740 + test_try("Test radius assingment after");
  741 + try {
  742 + GSkyRegionCircle circle(refdir_radeczerozero,refradius);
  743 + circle.radius(-1.);
  744 + test_try_failure();
  745 + }
  746 + catch (GException::invalid_argument &e) {
  747 + test_try_success();
  748 + }
  749 + catch (std::exception &e) {
  750 + test_try_failure(e);
  751 + }
  752 +
  753 + // Check radius assignment
  754 + GSkyRegionCircle refregion(refdir_radeczerozero,refradius);
  755 + double refradius_check = refregion.radius();
  756 + test_value(refradius,refradius_check,1.0e-10, "Test radius assignment");
  757 +
  758 + // Check solid angle assignment
  759 + double solidangle_check = refregion.solidangle();
  760 + double solidangle = 2*gammalib::pi*(1- std::cos(refradius /180 * gammalib::pi));
  761 + test_value(solidangle_check,solidangle,1.0e-10, "Test solid angle assignment");
  762 +
  763 + //exit test
  764 + return;
  765 +}
  766 +
  767 +/***************************************************************************
  768 + * @brief GSkyRegions_logic
  769 + ***************************************************************************/
  770 +void TestGSky::test_GSkyRegionCircle_logic(void)
  771 +{
  772 +
  773 + // set input files /filenames
  774 +// const std::string file1 = "test_GSkyRegionCircle1.reg";
  775 +// const std::string file2 = "test_GSkyRegionCircle2.reg";
  776 +
  777 + //set reference values
  778 +
  779 + GSkyDir refdir_radeczerozero = GSkyDir();
  780 + refdir_radeczerozero.radec_deg(0,0);
  781 +
  782 + GSkyDir refdir_lbzerozero = GSkyDir();
  783 + refdir_lbzerozero.lb_deg(0,0);
  784 +
  785 + GSkyDir refdir_outside_refregion = GSkyDir();
  786 + refdir_outside_refregion.radec_deg(100,80);
  787 +
  788 + GSkyDir refdir_raoffset = GSkyDir();
  789 + refdir_raoffset.radec_deg(5,0);
  790 +
  791 + GSkyDir refdir_decpole = GSkyDir();
  792 + refdir_decpole.radec_deg(0,89);
  793 +
  794 + GSkyDir refdir_ndecpole = GSkyDir();
  795 + refdir_ndecpole.radec_deg(100,89);
  796 +
  797 + GSkyDir refdir_rapole = GSkyDir();
  798 + refdir_rapole.radec_deg(359,0);
  799 +
  800 + GSkyRegionCircle refregion(refdir_radeczerozero,10);
  801 + GSkyDir centre = refregion.centre();
  802 +
  803 + GSkyRegionCircle refregion_smaller(refdir_radeczerozero,5);
  804 + GSkyRegionCircle refregion_larger(refdir_radeczerozero,20);
  805 + GSkyRegionCircle refregion_raoffset(refdir_raoffset,10);
  806 + GSkyRegionCircle refregion_rapole(refdir_rapole,3);
  807 + GSkyRegionCircle refregion_decpole(refdir_decpole,3);
  808 +
  809 + // Test contain dirs
  810 + test_assert(refregion.contains(refdir_radeczerozero),"test for containment");
  811 + test_assert(refregion.contains(centre), "test if center in circle");
  812 + test_assert(!refregion.contains(refdir_outside_refregion), "test2 for containment");
  813 +
  814 + // Test contain regions
  815 + test_assert(refregion.contains(refregion_smaller),"test for containment region");
  816 + test_assert(!refregion.contains(refregion_larger), "test for containment region2 ");
  817 + test_assert(refregion.contains(refregion), "test3 for containment region");
  818 +
  819 + test_assert(refregion.contains(refdir_rapole), "rapole for containment region");
  820 + test_assert(refregion_decpole.contains(refdir_ndecpole), "rapole for containment region");
  821 +
  822 + // Test overlaps
  823 + test_assert(refregion.overlaps(refregion_smaller),"test for overlap");
  824 + test_assert(refregion.overlaps(refregion_larger),"test2 for overlap");
  825 + test_assert(refregion.overlaps(refregion_raoffset),"test3 for overlap");
  826 + test_assert(!refregion.overlaps(refregion_decpole),"test4 for overlap");
  827 +
  828 + // Exit test
  829 + return;
  830 +}
  831 +
  832 +
  833 +/***************************************************************************
  834 + * @brief Test GSkyRegions input and output
  835 + ***************************************************************************/
  836 +void TestGSky::test_GSkyRegions_io(void)
  837 +{
  838 + // Set filenames
  839 + const std::string filereg_radeczerozero = "data/test_circle_radeczerozero.reg";
  840 +
  841 + // Test regions loading
  842 + test_try("Test regions loading");
  843 + try {
  844 + GSkyRegions regions;
  845 + regions.load(filereg_radeczerozero);
  846 + test_try_success();
  847 + }
  848 + catch (std::exception &e) {
  849 + test_try_failure(e);
  850 + }
  851 +
  852 +// //check if load does the proper thing
  853 +// GSkyRegions regs;
  854 +// GSkyRegion* reg;
  855 +// GSkyRegionCircle* circ;
  856 +// GSkyDir dir;
  857 +// double rad;
  858 +//
  859 +// regs.load(filereg_radeczerozero);
  860 +// reg = regs[0];
  861 +// if (reg->type() == "circle"){
  862 +// circ = (GSkyRegionCircle*)reg;
  863 +// dir = circ->centre();
  864 +// rad = circ->radius();
  865 +// }
  866 +// GSkyDir refdir_radeczerozero = GSkyDir();
  867 +// refdir_radeczerozero.radec_deg(0,0);
  868 +//
  869 +// test_value(rad,10,1.0e-10, "Test radius assignment");
  870 +//
  871 +//
  872 +
  873 + // Test regions saving
  874 + test_try("Test regions saving");
  875 + try {
  876 + GSkyRegions regions;
  877 + regions.save("region.ds9");
  878 + test_try_success();
  879 + }
  880 + catch (std::exception &e) {
  881 + test_try_failure(e);
  882 + }
  883 +
  884 + // Exit test
  885 + return;
  886 +}
  887 +
  888 +
  889 +/***************************************************************************
  890 + * @brief Main test function
704 891 ***************************************************************************/
705 892 int main(void)
706 893 {
... ...
test/test_GSky.hpp
  View file @7656e43
... ... @@ -40,6 +40,9 @@ class TestGSky : public GTestSuite
40 40 void test_GSkymap_healpix_io(void);
41 41 void test_GSkymap_wcs_construct(void);
42 42 void test_GSkymap_wcs_io(void);
  43 + void test_GSkyRegions_io(void);
  44 + void test_GSkyRegionCircle_construct(void);
  45 + void test_GSkyRegionCircle_logic(void);
43 46  
44 47 // Private methods
45 48 private:
... ...
test/test_GXspec.cpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * test_GXspec.hpp - Test Xspec module *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file test_GXspec.cpp
  23 + * @brief Implementation of unit tests for Xspec module
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "test_GXspec.hpp"
  32 +#include "GTools.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @brief Set parameters and tests
  37 + **************************************************************************/
  38 +void TestGXspec::set(void)
  39 +{
  40 + // Test name
  41 + name("GXspec");
  42 +
  43 + // Append tests
  44 + append(static_cast<pfunction>(&TestGXspec::test_GPha), "Test GPha class");
  45 + append(static_cast<pfunction>(&TestGXspec::test_GArf), "Test GArf class");
  46 + append(static_cast<pfunction>(&TestGXspec::test_GRmf), "Test GRmf class");
  47 +
  48 + // Return
  49 + return;
  50 +}
  51 +
  52 +
  53 +/***********************************************************************//**
  54 + * @brief Test GPha class
  55 + **************************************************************************/
  56 +void TestGXspec::test_GPha(void)
  57 +{
  58 + // Test void constructor
  59 + test_try("GPha void constructor");
  60 + try {
  61 + GPha pha;
  62 + test_try_success();
  63 + }
  64 + catch (std::exception &e) {
  65 + test_try_failure(e);
  66 + }
  67 +
  68 + // Test energy boundary constructor
  69 + test_try("GPha energy boundary constructor");
  70 + try {
  71 + GEbounds ebds(10, GEnergy(0.1, "TeV"), GEnergy(10.0, "TeV"));
  72 + GPha pha(ebds);
  73 + test_try_success();
  74 + }
  75 + catch (std::exception &e) {
  76 + test_try_failure(e);
  77 + }
  78 +
  79 + // Test filling and accessing
  80 + GEbounds ebds(9, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"), false);
  81 + GPha pha(ebds);
  82 + for (int i = 0; i < 12; i += 2) {
  83 + GEnergy energy(double(i), "TeV");
  84 + pha.fill(energy);
  85 + }
  86 + test_value(pha.counts(), 4.0);
  87 + test_value(pha.underflow(), 1.0);
  88 + test_value(pha.overflow(), 1.0);
  89 + test_value(pha.outflow(), 0.0);
  90 + for (int i = 0; i < 9; i += 2) {
  91 + test_value(pha.at(i), 0.0);
  92 + test_value(pha[i], 0.0);
  93 + }
  94 + for (int i = 1; i < 9; i += 2) {
  95 + test_value(pha.at(i), 1.0);
  96 + test_value(pha[i], 1.0);
  97 + }
  98 + pha[0] = 5.0;
  99 + pha[1] = 3.7;
  100 + test_value(pha[0], 5.0);
  101 + test_value(pha[1], 3.7);
  102 +
  103 + // Test saving and loading
  104 + pha.save("pha.fits", true);
  105 + test_assert(pha.filename() == "pha.fits",
  106 + "Unexpected filename \""+pha.filename()+"\".");
  107 + pha.load("pha.fits");
  108 + test_assert(pha.filename() == "pha.fits",
  109 + "Unexpected filename \""+pha.filename()+"\".");
  110 + test_value(pha[0], 5.0, 1.0e-6);
  111 + test_value(pha[1], 3.7, 1.0e-6);
  112 + test_value(pha.counts(), 11.7, 1.0e-6);
  113 + test_value(pha.underflow(), 0.0, 1.0e-6);
  114 + test_value(pha.overflow(), 0.0, 1.0e-6);
  115 + test_value(pha.outflow(), 0.0, 1.0e-6);
  116 + for (int i = 2; i < 9; i += 2) {
  117 + test_value(pha.at(i), 0.0);
  118 + test_value(pha[i], 0.0);
  119 + }
  120 + for (int i = 3; i < 9; i += 2) {
  121 + test_value(pha.at(i), 1.0);
  122 + test_value(pha[i], 1.0);
  123 + }
  124 +
  125 + // Test constructing
  126 + GPha pha2("pha.fits");
  127 + test_assert(pha2.filename() == "pha.fits",
  128 + "Unexpected filename \""+pha2.filename()+"\".");
  129 + test_value(pha2[0], 5.0, 1.0e-6);
  130 + test_value(pha2[1], 3.7, 1.0e-6);
  131 + test_value(pha2.counts(), 11.7, 1.0e-6);
  132 + test_value(pha2.underflow(), 0.0, 1.0e-6);
  133 + test_value(pha2.overflow(), 0.0, 1.0e-6);
  134 + test_value(pha2.outflow(), 0.0, 1.0e-6);
  135 + for (int i = 2; i < 9; i += 2) {
  136 + test_value(pha2.at(i), 0.0);
  137 + test_value(pha2[i], 0.0);
  138 + }
  139 + for (int i = 3; i < 9; i += 2) {
  140 + test_value(pha2.at(i), 1.0);
  141 + test_value(pha2[i], 1.0);
  142 + }
  143 + //std::cout << pha << std::endl;
  144 + //std::cout << pha2 << std::endl;
  145 +
  146 + // Return
  147 + return;
  148 +}
  149 +
  150 +
  151 +/***********************************************************************//**
  152 + * @brief Test GArf class
  153 + **************************************************************************/
  154 +void TestGXspec::test_GArf(void)
  155 +{
  156 + // Test void constructor
  157 + test_try("GArf void constructor");
  158 + try {
  159 + GArf arf;
  160 + test_try_success();
  161 + }
  162 + catch (std::exception &e) {
  163 + test_try_failure(e);
  164 + }
  165 +
  166 + // Test energy boundary constructor
  167 + test_try("GArf energy boundary constructor");
  168 + try {
  169 + GEbounds ebds(10, GEnergy(0.1, "TeV"), GEnergy(10.0, "TeV"));
  170 + GArf arf(ebds);
  171 + test_try_success();
  172 + }
  173 + catch (std::exception &e) {
  174 + test_try_failure(e);
  175 + }
  176 +
  177 + // Test filling and accessing
  178 + GEbounds ebds(9, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"));
  179 + GArf arf(ebds);
  180 + for (int i = 0; i < 9; i += 2) {
  181 + arf[i] = 1.0;
  182 + }
  183 + for (int i = 0; i < 9; i += 2) {
  184 + test_value(arf.at(i), 1.0);
  185 + test_value(arf[i], 1.0);
  186 + }
  187 + for (int i = 1; i < 9; i += 2) {
  188 + test_value(arf.at(i), 0.0);
  189 + test_value(arf[i], 0.0);
  190 + }
  191 + arf[0] = 5.0;
  192 + arf[1] = 3.7;
  193 + test_value(arf[0], 5.0);
  194 + test_value(arf[1], 3.7);
  195 +
  196 + // Test saving and loading
  197 + arf.save("arf.fits", true);
  198 + test_assert(arf.filename() == "arf.fits",
  199 + "Unexpected filename \""+arf.filename()+"\".");
  200 + arf.load("arf.fits");
  201 + test_assert(arf.filename() == "arf.fits",
  202 + "Unexpected filename \""+arf.filename()+"\".");
  203 + test_value(arf[0], 5.0, 1.0e-6);
  204 + test_value(arf[1], 3.7, 1.0e-6);
  205 + for (int i = 2; i < 9; i += 2) {
  206 + test_value(arf.at(i), 1.0);
  207 + test_value(arf[i], 1.0);
  208 + }
  209 + for (int i = 3; i < 9; i += 2) {
  210 + test_value(arf.at(i), 0.0);
  211 + test_value(arf[i], 0.0);
  212 + }
  213 +
  214 + // Test constructing
  215 + GArf arf2("arf.fits");
  216 + test_assert(arf2.filename() == "arf.fits",
  217 + "Unexpected filename \""+arf2.filename()+"\".");
  218 + test_value(arf2[0], 5.0, 1.0e-6);
  219 + test_value(arf2[1], 3.7, 1.0e-6);
  220 + for (int i = 2; i < 9; i += 2) {
  221 + test_value(arf2.at(i), 1.0);
  222 + test_value(arf2[i], 1.0);
  223 + }
  224 + for (int i = 3; i < 9; i += 2) {
  225 + test_value(arf2.at(i), 0.0);
  226 + test_value(arf2[i], 0.0);
  227 + }
  228 + //std::cout << arf << std::endl;
  229 + //std::cout << arf2 << std::endl;
  230 +
  231 + // Return
  232 + return;
  233 +}
  234 +
  235 +
  236 +/***********************************************************************//**
  237 + * @brief Test GRmf class
  238 + **************************************************************************/
  239 +void TestGXspec::test_GRmf(void)
  240 +{
  241 + // Test void constructor
  242 + test_try("GRmf void constructor");
  243 + try {
  244 + GRmf rmf;
  245 + test_try_success();
  246 + }
  247 + catch (std::exception &e) {
  248 + test_try_failure(e);
  249 + }
  250 +
  251 + // Test energy boundary constructor
  252 + test_try("GRmf energy boundary constructor");
  253 + try {
  254 + GEbounds ebds(10, GEnergy(0.1, "TeV"), GEnergy(10.0, "TeV"));
  255 + GRmf rmf(ebds, ebds);
  256 + test_try_success();
  257 + }
  258 + catch (std::exception &e) {
  259 + test_try_failure(e);
  260 + }
  261 +
  262 + // Test filling and accessing
  263 + GEbounds ebds(9, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"));
  264 + GRmf rmf(ebds, ebds);
  265 + for (int i = 0; i < 9; ++i) {
  266 + for (int k = i; k < i+3 && k < 9; ++k) {
  267 + rmf(i, k) = 1.0;
  268 + }
  269 + }
  270 + for (int i = 0; i < 9; ++i) {
  271 + int k = 0;
  272 + for (; k < i; ++k) {
  273 + test_value(rmf.at(i,k), 0.0);
  274 + test_value(rmf(i,k), 0.0);
  275 + }
  276 + for (; k < i+3 && k < 9; ++k) {
  277 + test_value(rmf.at(i,k), 1.0);
  278 + test_value(rmf(i,k), 1.0);
  279 + }
  280 + for (; k < 9; ++k) {
  281 + test_value(rmf.at(i,k), 0.0);
  282 + test_value(rmf(i,k), 0.0);
  283 + }
  284 + }
  285 +
  286 + // Test saving and loading
  287 + rmf.save("rmf.fits", true);
  288 + test_assert(rmf.filename() == "rmf.fits",
  289 + "Unexpected filename \""+rmf.filename()+"\".");
  290 + rmf.load("rmf.fits");
  291 + test_assert(rmf.filename() == "rmf.fits",
  292 + "Unexpected filename \""+rmf.filename()+"\".");
  293 + /*
  294 + for (int i = 0; i < 9; ++i) {
  295 + for (int k = 0; k < 9; ++k) {
  296 + std::cout << rmf(i,k) << " ";
  297 + }
  298 + std::cout << std::endl;
  299 + }
  300 + */
  301 + for (int i = 0; i < 9; ++i) {
  302 + int k = 0;
  303 + for (; k < i; ++k) {
  304 + test_value(rmf.at(i,k), 0.0);
  305 + test_value(rmf(i,k), 0.0);
  306 + }
  307 + for (; k < i+3 && k < 9; ++k) {
  308 + test_value(rmf.at(i,k), 1.0);
  309 + test_value(rmf(i,k), 1.0);
  310 + }
  311 + for (; k < 9; ++k) {
  312 + test_value(rmf.at(i,k), 0.0);
  313 + test_value(rmf(i,k), 0.0);
  314 + }
  315 + }
  316 +
  317 + // Test constructing
  318 + GRmf rmf2("rmf.fits");
  319 + test_assert(rmf2.filename() == "rmf.fits",
  320 + "Unexpected filename \""+rmf2.filename()+"\".");
  321 + for (int i = 0; i < 9; ++i) {
  322 + int k = 0;
  323 + for (; k < i; ++k) {
  324 + test_value(rmf.at(i,k), 0.0);
  325 + test_value(rmf(i,k), 0.0);
  326 + }
  327 + for (; k < i+3 && k < 9; ++k) {
  328 + test_value(rmf.at(i,k), 1.0);
  329 + test_value(rmf(i,k), 1.0);
  330 + }
  331 + for (; k < 9; ++k) {
  332 + test_value(rmf.at(i,k), 0.0);
  333 + test_value(rmf(i,k), 0.0);
  334 + }
  335 + }
  336 + //std::cout << rmf << std::endl;
  337 + //std::cout << rmf2 << std::endl;
  338 +
  339 + // Return
  340 + return;
  341 +}
  342 +
  343 +
  344 +/***************************************************************************
  345 + * @brief Main entry point for test executable
  346 + ***************************************************************************/
  347 +int main(void)
  348 +{
  349 + // Allocate test suit container
  350 + GTestSuites testsuites("Xspec module");
  351 +
  352 + // Initially assume that we pass all tests
  353 + bool success = true;
  354 +
  355 + // Create a test suite
  356 + TestGXspec test;
  357 +
  358 + // Append test suite to the container
  359 + testsuites.append(test);
  360 +
  361 + // Run the testsuites
  362 + success = testsuites.run();
  363 +
  364 + // Save test report
  365 + testsuites.save("reports/GXspec.xml");
  366 +
  367 + // Return success status
  368 + return (success ? 0 : 1);
  369 +}
... ...
test/test_GXspec.hpp 0 → 100644
  View file @7656e43
  1 +/***************************************************************************
  2 + * test_GXspec.hpp - Test Xspec module *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file test_GXspec.hpp
  23 + * @brief Definition of unit tests for Xspec module
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef TEST_GXSPEC_HPP
  28 +#define TEST_GXSPEC_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GammaLib.hpp"
  32 +
  33 +
  34 +/***********************************************************************//**
  35 + * @class TestGXspec
  36 + *
  37 + * @brief Test suite for Xspec module
  38 + ***************************************************************************/
  39 +class TestGXspec : public GTestSuite {
  40 +
  41 +public:
  42 + // Constructors and destructors
  43 + TestGXspec(void) : GTestSuite() {}
  44 + virtual ~TestGXspec(void) {}
  45 +
  46 + // Methods
  47 + virtual void set(void);
  48 + void test_GPha(void);
  49 + void test_GArf(void);
  50 + void test_GRmf(void);
  51 +
  52 +private:
  53 + // Private members
  54 +};
  55 +
  56 +#endif /* TEST_GXSPEC_HPP */
... ...
test/test_GXspec.py 0 → 100755
  View file @7656e43
  1 +# ==========================================================================
  2 +# This module performs unit tests for the GammaLib xspec module.
  3 +#
  4 +# Copyright (C) 2013 Juergen Knoedlseder
  5 +#
  6 +# This program is free software: you can redistribute it and/or modify
  7 +# it under the terms of the GNU General Public License as published by
  8 +# the Free Software Foundation, either version 3 of the License, or
  9 +# (at your option) any later version.
  10 +#
  11 +# This program is distributed in the hope that it will be useful,
  12 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  14 +# GNU General Public License for more details.
  15 +#
  16 +# You should have received a copy of the GNU General Public License
  17 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  18 +#
  19 +# ==========================================================================
  20 +from gammalib import *
  21 +from math import *
  22 +import os
  23 +
  24 +
  25 +# ==================================== #
  26 +# Test class for GammaLib xspec module #
  27 +# ==================================== #
  28 +class Test(GPythonTestSuite):
  29 + """
  30 + Test class for GammaLib xspec module.
  31 + """
  32 + # Constructor
  33 + def __init__(self):
  34 + """
  35 + Constructor.
  36 + """
  37 + # Call base class constructor
  38 + GPythonTestSuite.__init__(self)
  39 +
  40 + # Return
  41 + return
  42 +
  43 + # Set test functions
  44 + def set(self):
  45 + """
  46 + Set all test functions.
  47 + """
  48 + # Set test name
  49 + self.name("xspec")
  50 +
  51 + # Append tests
  52 + self.append(self.test_pha, "Test GPha")
  53 + self.append(self.test_arf, "Test GArf")
  54 + self.append(self.test_rmf, "Test GRmf")
  55 +
  56 + # Return
  57 + return
  58 +
  59 + # Test GPha class
  60 + def test_pha(self):
  61 + """
  62 + Test GPha class.
  63 + """
  64 + # Allocate GPha object and check its size
  65 + nmeasured = 10
  66 + emeasured = GEbounds(nmeasured, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"))
  67 + pha = GPha(emeasured)
  68 + self.test_value(pha.size(), nmeasured)
  69 +
  70 + # Set and retrieve a value
  71 + value = 3.7
  72 + pha[3] = value
  73 + self.test_value(pha[3], value)
  74 +
  75 + # Return
  76 + return
  77 +
  78 + # Test GArf class
  79 + def test_arf(self):
  80 + """
  81 + Test GArf class.
  82 + """
  83 + # Allocate GArf object and check its size
  84 + nmeasured = 10
  85 + emeasured = GEbounds(nmeasured, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"))
  86 + arf = GArf(emeasured)
  87 + self.test_value(arf.size(), nmeasured)
  88 +
  89 + # Set and retrieve a value
  90 + value = 3.7
  91 + arf[3] = value
  92 + self.test_value(arf[3], value)
  93 +
  94 + # Return
  95 + return
  96 +
  97 + # Test GRmf class
  98 + def test_rmf(self):
  99 + """
  100 + Test GRmf class.
  101 + """
  102 + # Allocate GRmf object and check its size
  103 + ntrue = 20
  104 + nmeasured = 10
  105 + etrue = GEbounds(ntrue, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"))
  106 + emeasured = GEbounds(nmeasured, GEnergy(1.0, "TeV"), GEnergy(10.0, "TeV"))
  107 + rmf = GRmf(etrue, emeasured)
  108 + self.test_value(rmf.size(), ntrue*nmeasured)
  109 + self.test_value(rmf.ntrue(), ntrue)
  110 + self.test_value(rmf.nmeasured(), nmeasured)
  111 +
  112 + # Set and retrieve a value
  113 + value = 3.7
  114 + rmf[3,5] = value
  115 + self.test_value(rmf[3,5], value)
  116 +
  117 + # Return
  118 + return
... ...
test/test_python.py
  View file @7656e43
... ... @@ -32,6 +32,7 @@ import test_GSky
32 32 import test_GSupport
33 33 import test_GTest
34 34 import test_GXml
  35 +import test_GXspec
35 36  
36 37 # Try importing MWL tests
37 38 try:
... ... @@ -73,17 +74,18 @@ if __name__ == &#39;__main__&#39;:
73 74 suites = GTestSuites("Python interface testing")
74 75  
75 76 # Allocate test suites and append them to the container
76   - suite1 = test_GApplication.Test()
77   - suite2 = test_GFits.Test()
78   - suite3 = test_GLinalg.Test()
79   - suite4 = test_GModel.Test()
80   - suite5 = test_GNumerics.Test()
81   - suite6 = test_GObservation.Test()
82   - suite7 = test_GOptimizer.Test()
83   - suite8 = test_GSky.Test()
84   - suite9 = test_GSupport.Test()
  77 + suite1 = test_GApplication.Test()
  78 + suite2 = test_GFits.Test()
  79 + suite3 = test_GLinalg.Test()
  80 + suite4 = test_GModel.Test()
  81 + suite5 = test_GNumerics.Test()
  82 + suite6 = test_GObservation.Test()
  83 + suite7 = test_GOptimizer.Test()
  84 + suite8 = test_GSky.Test()
  85 + suite9 = test_GSupport.Test()
85 86 suite10 = test_GTest.Test()
86 87 suite11 = test_GXml.Test()
  88 + suite12 = test_GXspec.Test()
87 89  
88 90 # Setup unit tests
89 91 suite1.set()
... ... @@ -97,6 +99,7 @@ if __name__ == &#39;__main__&#39;:
97 99 suite9.set()
98 100 suite10.set()
99 101 suite11.set()
  102 + suite12.set()
100 103  
101 104 # Append tests to container
102 105 suites.append(suite1)
... ... @@ -110,6 +113,7 @@ if __name__ == &#39;__main__&#39;:
110 113 suites.append(suite9)
111 114 suites.append(suite10)
112 115 suites.append(suite11)
  116 + suites.append(suite12)
113 117  
114 118 # Optionally handle MWL suite
115 119 if has_mwl:
... ...