Stefan Eschbach / gammalib

Commit 956995475d7c4f1de1538fc26d57514c06501196

Authored by Jürgen Knödlseder
2 parents 986c4d09 82c2d3e1

Merge branch '1996-ctselect-add-phase-info-GCTAEventList' into integration 

# Conflicts:
#	ChangeLog
#	inst/cta/include/GCTAEventList.hpp
#	inst/cta/pyext/GCTAEventList.i
Inline Side-by-side
Showing 18 changed files with 900 additions and 14 deletions
ChangeLog
  View file @9569954
... ... @@ -3,6 +3,7 @@
3 3 * Version 1.3.0 released
4 4 =========================
5 5  
  6 + Add GPhases class (#1996)
6 7 Add support for Monte Carlo identifier metadata (#1483)
7 8 Add GEventList::has_mc_id() and GEventAtom::mc_id() methods (#1483)
8 9 Add GEventList::has_phase() and GEventList::has_detxy() methods (#1982)
... ...
... ... @@ -12,6 +12,7 @@ made since the last release of GammaLib.
12 12 2. Public interface modifications
13 13 ---------------------------------
14 14 The following classes have been added:
  15 +- GPhases
15 16 - GWcsGLS
16 17 - GWcsSFL
17 18 - GModelSpectralSmoothBrokenPlaw
... ... @@ -139,6 +140,8 @@ GEbounds::set_log() method throws an exception if either the minimum or
139 140 maximum energy are zero. Furthermore, the methods clear the energy boundaries
140 141 if the requested number of intervals is non-zero (#1989).
141 142  
  143 +The GPhases class has been added to handle phase selections of events (#1996).
  144 +
142 145  
143 146 11. Optimization module
144 147 -----------------------
... ...
include/GPhases.hpp 0 → 100644
  View file @9569954
  1 +/***************************************************************************
  2 + * GPhases.hpp - Phase intervals class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GPhases.hpp
  23 + * @brief Phase intervals class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GPHASES_HPP
  28 +#define GPHASES_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GContainer.hpp"
  33 +
  34 +/* __ Forward declarations _______________________________________________ */
  35 +
  36 +
  37 +/***********************************************************************//**
  38 + * @class GPhases
  39 + *
  40 + * @brief Phase Intervals class
  41 + ***************************************************************************/
  42 +class GPhases : public GContainer {
  43 +
  44 +public:
  45 + // Constructors and destructors
  46 + GPhases(void);
  47 + GPhases(const GPhases& phases);
  48 + GPhases(const double& pmin, const double& pmax);
  49 + virtual ~GPhases(void);
  50 +
  51 + // Operators
  52 + GPhases& operator=(const GPhases& phases);
  53 +
  54 + // Methods
  55 + void clear(void);
  56 + GPhases* clone(void) const;
  57 + std::string classname(void) const;
  58 + int size(void) const;
  59 + bool is_empty(void) const;
  60 + void append(const double& pmin, const double& pmax);
  61 + void remove(const int& index);
  62 + void reserve(const int& num);
  63 + void extend(const GPhases& phases);
  64 + double pmin(const int& index) const;
  65 + double pmax(const int& index) const;
  66 + std::string print(const GChatter& chatter = NORMAL) const;
  67 +
  68 +protected:
  69 + // Protected methods
  70 + void init_members(void);
  71 + void copy_members(const GPhases& phases);
  72 + void free_members(void);
  73 + void insert_interval(const int& index, const double& pmin,
  74 + const double& pmax);
  75 +
  76 + // Protected data area
  77 + std::vector<double> m_pmin; // Lower phase interval boundaries
  78 + std::vector<double> m_pmax; // Upper phase interval boundaries
  79 +};
  80 +
  81 +
  82 +/***********************************************************************//**
  83 + * @brief Return class name
  84 + *
  85 + * @return String containing the class name ("GPhases").
  86 + ***************************************************************************/
  87 +inline
  88 +std::string GPhases::classname(void) const
  89 +{
  90 + return ("GPhases");
  91 +}
  92 +
  93 +
  94 +/***********************************************************************//**
  95 + * @brief Return number of phase intervals
  96 + *
  97 + * @return Number of phase intervals.
  98 + ***************************************************************************/
  99 +inline
  100 +int GPhases::size(void) const
  101 +{
  102 + return (int)m_pmin.size();
  103 +}
  104 +
  105 +
  106 +/***********************************************************************//**
  107 + * @brief Signal if there are no phase intervals
  108 + *
  109 + * @return True if there are no phase intervals.
  110 + ***************************************************************************/
  111 +inline
  112 +bool GPhases::is_empty(void) const
  113 +{
  114 + return (m_pmin.size() == 0);
  115 +}
  116 +
  117 +#endif /* GPHASES_HPP */
... ...
include/GammaLib.hpp
  View file @9569954
... ... @@ -183,6 +183,7 @@
183 183 #include "GGti.hpp"
184 184 #include "GRoi.hpp"
185 185 #include "GEbounds.hpp"
  186 +#include "GPhases.hpp"
186 187 #include "GInstDir.hpp"
187 188 #include "GResponse.hpp"
188 189 #include "GPhotons.hpp"
... ...
include/Makefile.am
  View file @9569954
... ... @@ -143,6 +143,7 @@ pkginclude_HEADERS = GammaLib.hpp \
143 143 GGti.hpp \
144 144 GRoi.hpp \
145 145 GEbounds.hpp \
  146 + GPhases.hpp \
146 147 GInstDir.hpp \
147 148 GResponse.hpp \
148 149 GPhotons.hpp \
... ...
inst/cta/include/GCTAEventList.hpp
  View file @9569954
... ... @@ -31,6 +31,7 @@
31 31 #include <string>
32 32 #include <vector>
33 33 #include "GEventList.hpp"
  34 +#include "GPhases.hpp"
34 35 #include "GFilename.hpp"
35 36 #include "GCTAEventAtom.hpp"
36 37 #include "GCTARoi.hpp"
... ... @@ -96,6 +97,8 @@ public:
96 97 double irf_cache(const std::string& name, const int& index) const;
97 98 void irf_cache(const std::string& name, const int& index,
98 99 const double& irf) const;
  100 + const GPhases& phases(void) const;
  101 + void phases(const GPhases& phases);
99 102 const std::string& gtiname(void) const;
100 103 void has_phase(const bool& has_phase);
101 104 void has_detxy(const bool& has_detxy);
... ... @@ -124,6 +127,7 @@ protected:
124 127  
125 128 // Event list meta data
126 129 GCTARoi m_roi; //!< Region of interest
  130 + GPhases m_phases; //!< Phase intervals
127 131 int m_num_events; //!< Number of events
128 132 std::string m_gti_extname; //!< GTI extension name
129 133 bool m_has_phase; //!< Signal presence of phase
... ... @@ -192,6 +196,31 @@ const GCTARoi&amp; GCTAEventList::roi(void) const
192 196  
193 197  
194 198 /***********************************************************************//**
  199 + * @brief Return phase intervals
  200 + *
  201 + * @return Phase intervals.
  202 + ***************************************************************************/
  203 +inline
  204 +const GPhases& GCTAEventList::phases(void) const
  205 +{
  206 + return (m_phases);
  207 +}
  208 +
  209 +
  210 +/***********************************************************************//**
  211 + * @brief Set phase intervals
  212 + *
  213 + * @param[in] phases Phase intervals.
  214 + ***************************************************************************/
  215 +inline
  216 +void GCTAEventList::phases(const GPhases& phases)
  217 +{
  218 + m_phases = phases;
  219 + return;
  220 +}
  221 +
  222 +
  223 +/***********************************************************************//**
195 224 * @brief Reserves space for events
196 225 *
197 226 * @param[in] number Number of events.
... ...
inst/cta/pyext/GCTAEventList.i
  View file @9569954
... ... @@ -62,17 +62,21 @@ public:
62 62 void append_column(const GFitsTableCol& column);
63 63 void reserve(const int& number);
64 64 void remove(const int& index, const int& number = 1);
65   - void write(GFits& fits, const std::string& evtname,
66   - const std::string& gtiname) const;
67   - void fetch(void);
68   - void dispose(void);
69   - double irf_cache(const std::string& name, const int& index) const;
70   - void irf_cache(const std::string& name, const int& index,
  65 + void write(GFits& fits,
  66 + const std::string& evtname,
  67 + const std::string& gtiname) const;
  68 + void fetch(void) const;
  69 + void dispose(void) const;
  70 + double irf_cache(const std::string& name,
  71 + const int& index) const;
  72 + void irf_cache(const std::string& name,
  73 + const int& index,
71 74 const double& irf) const;
  75 + const GPhases& phases(void) const;
  76 + void phases(const GPhases& phases);
72 77 const std::string& gtiname(void) const;
73 78 void has_phase(const bool& has_phase);
74 79 void has_detxy(const bool& has_detxy);
75   - void has_mc_id(const bool& has_mc_id);
76 80 const bool& has_phase() const;
77 81 const bool& has_detxy() const;
78 82 const bool& has_mc_id() const;
... ...
inst/cta/src/GCTAEventList.cpp
  View file @9569954
... ... @@ -53,6 +53,8 @@
53 53 #define G_SET_MC_ID_NAMES "GCTAEventList::set_mc_id_names(std::vector<int>&,"\
54 54 " std::vector<std::string>&)"
55 55 #define G_FETCH "GCTAEventList::fetch()"
  56 +#define G_WRITE_DS_KEYS "GCTAEventList::write_ds_keys(GFitsHDU&, "\
  57 + "const std::string&) const"
56 58  
57 59 /* __ Macros _____________________________________________________________ */
58 60  
... ... @@ -439,6 +441,9 @@ void GCTAEventList::read(const GFits&amp; fits)
439 441 // Read energy boundaries from data sub-space keyword
440 442 m_ebounds = gammalib::read_ds_ebounds(events);
441 443  
  444 + // Read phase intervals from data sub-space keyword
  445 + m_phases = gammalib::read_ds_phase(events);
  446 +
442 447 // If the file name contains an expression then load the events now.
443 448 // We do this at the end to make sure that the GTI has been set before.
444 449 if (m_filename.has_expression()) {
... ... @@ -957,6 +962,7 @@ void GCTAEventList::init_members(void)
957 962 {
958 963 // Initialise members
959 964 m_roi.clear();
  965 + m_phases.clear();
960 966 m_events.clear();
961 967 m_columns.clear();
962 968 m_filename.clear();
... ... @@ -986,6 +992,7 @@ void GCTAEventList::copy_members(const GCTAEventList&amp; list)
986 992 {
987 993 // Copy members
988 994 m_roi = list.m_roi;
  995 + m_phases = list.m_phases;
989 996 m_events = list.m_events;
990 997 m_filename = list.m_filename;
991 998 m_num_events = list.m_num_events;
... ... @@ -1337,6 +1344,9 @@ void GCTAEventList::write_events(GFitsBinTable&amp; hdu) const
1337 1344 * @param[in] hdu FITS HDU.
1338 1345 * @param[in] gtiname Good Time Interval FITS extension.
1339 1346 *
  1347 + * @exception GException::vector_mismatch
  1348 + * Lower and upper phase boundaries have different size.
  1349 + *
1340 1350 * Writes the data sub-space keywords for an event list into the FITS HDU.
1341 1351 * The following keywords will be written:
1342 1352 *
... ... @@ -1410,6 +1420,30 @@ void GCTAEventList::write_ds_keys(GFitsHDU&amp; hdu, const std::string&amp; gtiname) con
1410 1420  
1411 1421 } // endif: RoI was valid
1412 1422  
  1423 + // Check if there are phase interval cuts
  1424 + if (m_phases.size() > 0) {
  1425 +
  1426 + // Set phase interval string
  1427 + std::string dsval4 = "";
  1428 + for (int i = 0; i < m_phases.size(); ++i) {
  1429 + if (i != 0) {
  1430 + dsval4 += ",";
  1431 + }
  1432 + dsval4 += gammalib::str(m_phases.pmin(i)) + ":" +
  1433 + gammalib::str(m_phases.pmax(i));
  1434 + }
  1435 +
  1436 + // Write DS keywords
  1437 + hdu.card("DSTYP"+gammalib::str(ndskeys), "PHASE",
  1438 + "Data sub-space type");
  1439 + hdu.card("DSUNI"+gammalib::str(ndskeys), "DIMENSIONLESS",
  1440 + "Data sub-space unit");
  1441 + hdu.card("DSVAL"+gammalib::str(ndskeys), dsval4,
  1442 + "Data sub-space value");
  1443 + ndskeys++;
  1444 +
  1445 + } // endif: there were phase interval cuts
  1446 +
1413 1447 // Set number of data selection keys
1414 1448 hdu.card("NDSKEYS", ndskeys, "Number of data sub-space keys");
1415 1449  
... ...
inst/cta/src/GCTASupport.cpp
  View file @9569954
1 1 /***************************************************************************
2 2 * GCTASupport.cpp - CTA support functions *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -36,13 +36,15 @@
36 36 #include "GMath.hpp"
37 37 #include "GFitsHDU.hpp"
38 38 #include "GEbounds.hpp"
  39 +#include "GPhases.hpp"
39 40 #include "GCTARoi.hpp"
40 41 #include "GCTAInstDir.hpp"
41 42  
42 43 /* __ Method name definitions ____________________________________________ */
43 44 #define G_READ_DS_ROI "gammalib::read_ds_roi(GFitsHDU&)"
44 45 #define G_READ_DS_EBOUNDS "gammalib::read_ds_ebounds(GFitsHDU&)"
45   -#define G_READ_DS_GTI "GGti gammalib::read_ds_gti_extname(GFitsHDU&)"
  46 +#define G_READ_DS_PHASE "gammalib::read_ds_phase(GFitsHDU&)"
  47 +#define G_READ_DS_GTI "gammalib::read_ds_gti_extname(GFitsHDU&)"
46 48  
47 49 /* __ Coding definitions _________________________________________________ */
48 50  
... ... @@ -259,6 +261,69 @@ GEbounds gammalib::read_ds_ebounds(const GFitsHDU&amp; hdu)
259 261  
260 262  
261 263 /***********************************************************************//**
  264 + * @brief Read phase boundary data sub-space keywords
  265 + *
  266 + * @param[in] hdu FITS HDU
  267 + * @return Phase intervals
  268 + *
  269 + * @exception GException::invalid_value
  270 + * Invalid phase data sub-space encountered
  271 + *
  272 + * Reads the phase boundary data sub-space keywords by searching for a
  273 + * DSTYPx keyword named "PHASE". The data sub-space information is expected
  274 + * to be in the format "0.1:0.3,0.5:0.7", where each subset of numbers
  275 + * represents the minimum and maximum phase. No detailed syntax
  276 + * checking is performed.
  277 + ***************************************************************************/
  278 +GPhases gammalib::read_ds_phase(const GFitsHDU& hdu)
  279 +{
  280 + // Initialise phase intervals
  281 + GPhases phases;
  282 +
  283 + // Get number of data sub-space keywords (default to 0 if keyword is
  284 + // not found)
  285 + int ndskeys = (hdu.has_card("NDSKEYS")) ? hdu.integer("NDSKEYS") : 0;
  286 +
  287 + // Loop over all data sub-space keys
  288 + for (int i = 1; i <= ndskeys; ++i) {
  289 +
  290 + // Set data sub-space key strings
  291 + std::string type_key = "DSTYP"+gammalib::str(i);
  292 + std::string unit_key = "DSUNI"+gammalib::str(i);
  293 + std::string value_key = "DSVAL"+gammalib::str(i);
  294 +
  295 + // Continue only if type_key is found and if this key is PHASE
  296 + if (hdu.has_card(type_key) && hdu.string(type_key) == "PHASE") {
  297 +
  298 + // Extract phase boundaries
  299 + std::string value = hdu.string(value_key);
  300 + std::vector<std::string> intervals = gammalib::split(value, ",");
  301 + for (int j = 0; j < intervals.size(); ++j) {
  302 + std::vector<std::string> values = gammalib::split(intervals[j], ":");
  303 + if (values.size() == 2) {
  304 + double phmin = gammalib::todouble(values[0]);
  305 + double phmax = gammalib::todouble(values[1]);
  306 + phases.append(phmin, phmax);
  307 + }
  308 + else {
  309 + std::string msg = "Invalid phase value \""+value+
  310 + "\" encountered in data selection key \""+
  311 + value_key+"\"";
  312 + throw GException::invalid_value(G_READ_DS_PHASE, msg);
  313 + }
  314 +
  315 + } //endfor: looped over phase intervals
  316 +
  317 + } // endif: PHASE type_key found
  318 +
  319 + } // endfor: looped over data selection keys
  320 +
  321 + // Return
  322 + return phases;
  323 +}
  324 +
  325 +
  326 +/***********************************************************************//**
262 327 * @brief Return Good Time Intervals extension name from data sub-space keywords
263 328 *
264 329 * @param[in] hdu FITS HDU
... ...
inst/cta/src/GCTASupport.hpp
  View file @9569954
1 1 /***************************************************************************
2 2 * GCTASupport.hpp - CTA support functions *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -29,6 +29,7 @@
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
  32 +#include <vector>
32 33  
33 34 /* __ Namespaces _________________________________________________________ */
34 35  
... ... @@ -38,6 +39,7 @@
38 39 class GFitsHDU;
39 40 class GCTARoi;
40 41 class GEbounds;
  42 +class GPhases;
41 43  
42 44 /* __ Prototypes _________________________________________________________ */
43 45 namespace gammalib {
... ... @@ -46,6 +48,7 @@ namespace gammalib {
46 48 const double& roi, const double& cosroi);
47 49 GCTARoi read_ds_roi(const GFitsHDU& hdu);
48 50 GEbounds read_ds_ebounds(const GFitsHDU& hdu);
  51 + GPhases read_ds_phase(const GFitsHDU& hdu);
49 52 std::string read_ds_gti_extname(const GFitsHDU& hdu);
50 53 }
51 54  
... ...
pyext/GGti.i
  View file @9569954
... ... @@ -35,7 +35,7 @@
35 35 *
36 36 * @brief Good Time Interval class
37 37 ***************************************************************************/
38   -class GGti : public GBase {
  38 +class GGti : public GContainer {
39 39  
40 40 public:
41 41 // Constructors and destructors
... ...
pyext/GPhases.i 0 → 100644
  View file @9569954
  1 +/***************************************************************************
  2 + * GPhases.i - Phase intervals class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GPhases.i
  23 + * @brief Phase intervals class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GPhases.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 + * @class GPhases
  35 + *
  36 + * @brief Phase intervals class
  37 + ***************************************************************************/
  38 +class GPhases : public GContainer {
  39 +
  40 +public:
  41 + // Constructors and destructors
  42 + GPhases(void);
  43 + GPhases(const GPhases& phases);
  44 + GPhases(const double& pmin, const double& pmax);
  45 + virtual ~GPhases(void);
  46 +
  47 + // Methods
  48 + void clear(void);
  49 + GPhases* clone(void) const;
  50 + std::string classname(void) const;
  51 + int size(void) const;
  52 + bool is_empty(void) const;
  53 + void append(const double& pmin, const double& pmax);
  54 + void remove(const int& index);
  55 + void reserve(const int& num);
  56 + void extend(const GPhases& phases);
  57 + double pmin(const int& index) const;
  58 + double pmax(const int& index) const;
  59 +};
  60 +
  61 +
  62 +/***********************************************************************//**
  63 + * @brief GPhases class extension
  64 + ***************************************************************************/
  65 +%extend GPhases {
  66 + int __len__() {
  67 + return (self->size());
  68 + }
  69 + GPhases copy() {
  70 + return (*self);
  71 + }
  72 +};
... ...
pyext/gammalib/obs.i
  View file @9569954
1 1 /***************************************************************************
2 2 * obs.i - Observation module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -171,6 +171,7 @@
171 171 %include "GEventBin.i"
172 172 %include "GInstDir.i"
173 173 %include "GRoi.i"
  174 +%include "GPhases.i"
174 175 %include "GResponse.i"
175 176 %include "GPhoton.i"
176 177 %include "GPhotons.i"
... ...
src/obs/GGti.cpp
  View file @9569954
... ... @@ -100,7 +100,7 @@ GGti::GGti(const GFilename&amp; filename)
100 100 *
101 101 * @param[in] gti Good Time Intervals.
102 102 *
103   - * Constructs Good Time Intervals by coping other Good Time Intervals.
  103 + * Constructs Good Time Intervals by copying other Good Time Intervals.
104 104 ***************************************************************************/
105 105 GGti::GGti(const GGti& gti)
106 106 {
... ...
src/obs/GPhases.cpp 0 → 100644
  View file @9569954
  1 +/***************************************************************************
  2 + * GPhases.cpp - Phase intervals class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 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 GPhases.cpp
  23 + * @brief Phase intervals class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GTools.hpp"
  33 +#include "GPhases.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_REMOVE "GPhases::remove(int&)"
  37 +#define G_PMIN "GPhases::pmin(int&)"
  38 +#define G_PMAX "GPhases::pmax(int&)"
  39 +#define G_INSERT_INTERVAL "GPhases::insert_interval(int&, double&, double&)"
  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 + * Constructs empty phase intervals.
  58 + ***************************************************************************/
  59 +GPhases::GPhases(void)
  60 +{
  61 + // Initialise class members
  62 + init_members();
  63 +
  64 + // Return
  65 + return;
  66 +}
  67 +
  68 +
  69 +/***********************************************************************//**
  70 + * @brief Copy constructor
  71 + *
  72 + * @param[in] phases Phase intervals.
  73 + *
  74 + * Constructs phase intervals by copying other phase intervals.
  75 + ***************************************************************************/
  76 +GPhases::GPhases(const GPhases& phases)
  77 +{
  78 + // Initialise class members
  79 + init_members();
  80 +
  81 + // Copy members
  82 + copy_members(phases);
  83 +
  84 + // Return
  85 + return;
  86 +}
  87 +
  88 +
  89 +/***********************************************************************//**
  90 + * @brief Single phase interval constructor
  91 + *
  92 + * @param[in] pmin Lower boundary of phase interval.
  93 + * @param[in] pmax Upper boundary of phase interval.
  94 + *
  95 + * Constructs Good Time Intervals from a single phase interval, given by
  96 + * [p@ pmin, @p pmax].
  97 + ***************************************************************************/
  98 +GPhases::GPhases(const double& pmin, const double& pmax)
  99 +{
  100 + // Initialise members
  101 + init_members();
  102 +
  103 + // Append phase interval
  104 + append(pmin, pmax);
  105 +
  106 + // Return
  107 + return;
  108 +}
  109 +
  110 +
  111 +/***********************************************************************//**
  112 + * @brief Destructor
  113 + ***************************************************************************/
  114 +GPhases::~GPhases(void)
  115 +{
  116 + // Free members
  117 + free_members();
  118 +
  119 + // Return
  120 + return;
  121 +}
  122 +
  123 +
  124 +/*==========================================================================
  125 + = =
  126 + = Operators =
  127 + = =
  128 + ==========================================================================*/
  129 +
  130 +/***********************************************************************//**
  131 + * @brief Assignment operator
  132 + *
  133 + * @param[in] phases Phase intervals.
  134 + * @return Phase intervals.
  135 + ***************************************************************************/
  136 +GPhases& GPhases::operator=(const GPhases& phases)
  137 +{
  138 + // Execute only if object is not identical
  139 + if (this != &phases) {
  140 +
  141 + // Free members
  142 + free_members();
  143 +
  144 + // Initialise private members
  145 + init_members();
  146 +
  147 + // Copy members
  148 + copy_members(phases);
  149 +
  150 + } // endif: object was not identical
  151 +
  152 + // Return this object
  153 + return *this;
  154 +}
  155 +
  156 +
  157 +/*==========================================================================
  158 + = =
  159 + = Public methods =
  160 + = =
  161 + ==========================================================================*/
  162 +
  163 +/***********************************************************************//**
  164 + * @brief Clear phase intervals
  165 + ***************************************************************************/
  166 +void GPhases::clear(void)
  167 +{
  168 + // Free members
  169 + free_members();
  170 +
  171 + // Initialise private members
  172 + init_members();
  173 +
  174 + // Return
  175 + return;
  176 +}
  177 +
  178 +
  179 +/***********************************************************************//**
  180 + * @brief Clone phase intervals
  181 + *
  182 + * @return Pointer to deep copy of phase intervals.
  183 + ***************************************************************************/
  184 +GPhases* GPhases::clone(void) const
  185 +{
  186 + return new GPhases(*this);
  187 +}
  188 +
  189 +
  190 +/***********************************************************************//**
  191 + * @brief Append phase interval
  192 + *
  193 + * @param[in] pmin Lower boundary of phase interval.
  194 + * @param[in] pmax Upper boundary of phase interval.
  195 + *
  196 + * Appends a phase interval at the end of the container.
  197 + ***************************************************************************/
  198 +void GPhases::append(const double& pmin, const double& pmax)
  199 +{
  200 + // Insert phase interval at the end of the list
  201 + insert_interval(size(), pmin, pmax);
  202 +
  203 + // Return
  204 + return;
  205 +}
  206 +
  207 +
  208 +/***********************************************************************//**
  209 + * @brief Remove phase interval
  210 + *
  211 + * @param[in] index Phase interval index (0,...,size()-1).
  212 + *
  213 + * Removes phase interval at @p index from the container. All intervals
  214 + * after the specified @p index are moved forward by one position.
  215 + ***************************************************************************/
  216 +void GPhases::remove(const int& index)
  217 +{
  218 + #if defined(G_RANGE_CHECK)
  219 + // If index is outside boundary then throw an error
  220 + if (index < 0 || index >= size()) {
  221 + throw GException::out_of_range(G_REMOVE, "Phase interval", index,
  222 + size());
  223 + }
  224 + #endif
  225 +
  226 + // Remove lower and upper boundaries
  227 + m_pmin.erase(m_pmin.begin() + index);
  228 + m_pmax.erase(m_pmax.begin() + index);
  229 +
  230 + // Return
  231 + return;
  232 +}
  233 +
  234 +
  235 +/***********************************************************************//**
  236 + * @brief Reserve space for phase intervals
  237 + *
  238 + * @param[in] num Number of elements.
  239 + ***************************************************************************/
  240 +void GPhases::reserve(const int& num)
  241 +{
  242 + // Reserves memory
  243 + m_pmin.reserve(num);
  244 + m_pmax.reserve(num);
  245 +
  246 + // Return
  247 + return;
  248 +}
  249 +
  250 +
  251 +/***********************************************************************//**
  252 + * @brief Append phase intervals
  253 + *
  254 + * @param[in] phases Phase intervals.
  255 + *
  256 + * Append phase intervals to the container.
  257 + ***************************************************************************/
  258 +void GPhases::extend(const GPhases& phases)
  259 +{
  260 + // Continue only if phase intervals are not empty
  261 + if (!phases.is_empty()) {
  262 +
  263 + // Get size. Note that we extract the size first to avoid an
  264 + // endless loop that arises when a container is appended to
  265 + // itself.
  266 + int num = phases.size();
  267 +
  268 + // Resize the boundary vectors
  269 + reserve(size() + num);
  270 +
  271 + // Loop over all elements and append them to container
  272 + for (int i = 0; i < num; ++i) {
  273 + m_pmin.push_back(phases.pmin(i));
  274 + m_pmax.push_back(phases.pmax(i));
  275 + }
  276 +
  277 + } // endif: Phase intervals were not empty
  278 +
  279 + // Return
  280 + return;
  281 +}
  282 +
  283 +
  284 +/***********************************************************************//**
  285 + * @brief Returns lower boundary for a given phase interval
  286 + *
  287 + * @param[in] index Phase interval index (0,...,size()-1).
  288 + * @return Lower boundary of phase interval.
  289 + *
  290 + * @exception GException::out_of_range
  291 + * Specified index is out of range.
  292 + ***************************************************************************/
  293 +double GPhases::pmin(const int& index) const
  294 +{
  295 + #if defined(G_RANGE_CHECK)
  296 + // If index is outside boundary then throw an error
  297 + if (index < 0 || index >= size()) {
  298 + throw GException::out_of_range(G_PMIN, "Phase interval", index,
  299 + size());
  300 + }
  301 + #endif
  302 +
  303 + // Return
  304 + return (m_pmin[index]);
  305 +}
  306 +
  307 +
  308 +/***********************************************************************//**
  309 + * @brief Returns upper boundary for a given phase interval
  310 + *
  311 + * @param[in] index Phase interval index (0,...,size()-1).
  312 + * @return Upper boundary of phase interval.
  313 + *
  314 + * @exception GException::out_of_range
  315 + * Specified index is out of range.
  316 + ***************************************************************************/
  317 +double GPhases::pmax(const int& index) const
  318 +{
  319 + #if defined(G_RANGE_CHECK)
  320 + // If index is outside boundary then throw an error
  321 + if (index < 0 || index >= size()) {
  322 + throw GException::out_of_range(G_PMAX, "Phase interval", index,
  323 + size());
  324 + }
  325 + #endif
  326 +
  327 + // Return
  328 + return (m_pmax[index]);
  329 +}
  330 +
  331 +
  332 +/***********************************************************************//**
  333 + * @brief Print phase intervals
  334 + *
  335 + * @param[in] chatter Chattiness.
  336 + * @return String containing phase interval information.
  337 + ***************************************************************************/
  338 +std::string GPhases::print(const GChatter& chatter) const
  339 +{
  340 + // Initialise result string
  341 + std::string result;
  342 +
  343 + // Continue only if chatter is not silent
  344 + if (chatter != SILENT) {
  345 +
  346 + // Append header
  347 + result.append("=== GPhases ===");
  348 +
  349 + // Append phase interval information
  350 + result.append("\n"+gammalib::parformat("Number of intervals"));
  351 + result.append(gammalib::str(size()));
  352 +
  353 + } // endif: chatter was not silent
  354 +
  355 + // Return result
  356 + return result;
  357 +}
  358 +
  359 +
  360 +/*==========================================================================
  361 + = =
  362 + = Private methods =
  363 + = =
  364 + ==========================================================================*/
  365 +
  366 +/***********************************************************************//**
  367 + * @brief Initialise class members
  368 + ***************************************************************************/
  369 +void GPhases::init_members(void)
  370 +{
  371 + // Initialise members
  372 + m_pmin.clear();
  373 + m_pmax.clear();
  374 +
  375 + // Return
  376 + return;
  377 +}
  378 +
  379 +
  380 +/***********************************************************************//**
  381 + * @brief Copy class members
  382 + *
  383 + * @param[in] phases Phase intervals.
  384 + ***************************************************************************/
  385 +void GPhases::copy_members(const GPhases& phases)
  386 +{
  387 + // Copy attributes
  388 + m_pmin = phases.m_pmin;
  389 + m_pmax = phases.m_pmax;
  390 +
  391 + // Return
  392 + return;
  393 +}
  394 +
  395 +
  396 +/***********************************************************************//**
  397 + * @brief Delete class members
  398 + ***************************************************************************/
  399 +void GPhases::free_members(void)
  400 +{
  401 + // Return
  402 + return;
  403 +}
  404 +
  405 +
  406 +/***********************************************************************//**
  407 + * @brief Insert phase interval
  408 + *
  409 + * @param[in] index Index after which interval is inserted.
  410 + * @param[in] pmin Lower boundary of interval.
  411 + * @param[in] pmax Upper boundary of interval.
  412 + *
  413 + * @exception GException::invalid_argument
  414 + * Invalid phase interval boundaries specified
  415 + *
  416 + * Inserts a phase interval before the interval with the specified @p index.
  417 + * If no interval with the specified index exists then append the interval
  418 + * at the end of the existing phase intervals.
  419 + ***************************************************************************/
  420 +void GPhases::insert_interval(const int& index,
  421 + const double& pmin,
  422 + const double& pmax)
  423 +{
  424 + // Throw an exception if phase minimum is outside the range [0,1]
  425 + if (pmin < 0.0 || pmin > 1.0) {
  426 + std::string msg = "Phase minimum "+gammalib::str(pmin)+
  427 + " outside the valid range [0,1]. Please "
  428 + "specify phase interval boundaries comprised "
  429 + "within 0 and 1.";
  430 + throw GException::invalid_argument(G_INSERT_INTERVAL, msg);
  431 + }
  432 +
  433 + // Throw an exception if phase minimum is outside the range [0,1]
  434 + if (pmax < 0.0 || pmax > 1.0) {
  435 + std::string msg = "Phase maximum "+gammalib::str(pmax)+
  436 + " outside the valid range [0,1]. Please "
  437 + "specify phase interval boundaries comprised "
  438 + "within 0 and 1.";
  439 + throw GException::invalid_argument(G_INSERT_INTERVAL, msg);
  440 + }
  441 +
  442 + // If phase minimum is smaller than phase maximum then append or insert
  443 + // a single interval ...
  444 + if (pmin < pmax) {
  445 + if (index >= size()) {
  446 + m_pmin.push_back(pmin);
  447 + m_pmax.push_back(pmax);
  448 + }
  449 + else {
  450 + m_pmin.insert(m_pmin.begin()+index, pmin);
  451 + m_pmax.insert(m_pmax.begin()+index, pmax);
  452 + }
  453 + }
  454 +
  455 + // ... otherwise if the phase minimum is larger than the maximum then
  456 + // consider this a wrap around interval and append or insert one interval
  457 + // from phase minimum to 1 and another interval from 0 to phase maximum
  458 + else if (pmin > pmax) {
  459 + if (pmin < 1.0) {
  460 + if (index >= size()) {
  461 + m_pmin.push_back(pmin);
  462 + m_pmax.push_back(1.0);
  463 + }
  464 + else {
  465 + m_pmin.insert(m_pmin.begin()+index, pmin);
  466 + m_pmax.insert(m_pmax.begin()+index, 1.0);
  467 + }
  468 + }
  469 + if (pmax > 0.0) {
  470 + if (index >= size()) {
  471 + m_pmin.push_back(0.0);
  472 + m_pmax.push_back(pmax);
  473 + }
  474 + else {
  475 + m_pmin.insert(m_pmin.begin()+index, 0.0);
  476 + m_pmax.insert(m_pmax.begin()+index, pmax);
  477 + }
  478 + }
  479 + }
  480 +
  481 + // Return
  482 + return;
  483 +}
... ...
src/obs/Makefile.am
  View file @9569954
1 1 #############################################################################
2 2 # Makefile.am for obs module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2017 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -44,6 +44,7 @@ sources = GEnergies.cpp \
44 44 GGti.cpp \
45 45 GRoi.cpp \
46 46 GEbounds.cpp \
  47 + GPhases.cpp \
47 48 GResponse.cpp \
48 49 GInstDir.cpp \
49 50 GPhotons.cpp \
... ...
test/test_GObservation.cpp
  View file @9569954
... ... @@ -69,6 +69,8 @@ void TestGObservation::set(void)
69 69 "Test GEnergies class");
70 70 append(static_cast<pfunction>(&TestGObservation::test_ebounds),
71 71 "Test GEbounds class");
  72 + append(static_cast<pfunction>(&TestGObservation::test_phases),
  73 + "Test GPhases class");
72 74 append(static_cast<pfunction>(&TestGObservation::test_photons),
73 75 "Test GPhotons class");
74 76 append(static_cast<pfunction>(&TestGObservation::test_observations_optimizer),
... ... @@ -1343,6 +1345,74 @@ void TestGObservation::test_photons(void)
1343 1345  
1344 1346  
1345 1347 /***********************************************************************//**
  1348 + * @brief Test GPhases
  1349 + ***************************************************************************/
  1350 +void TestGObservation::test_phases(void)
  1351 +{
  1352 + // Test void constructor
  1353 + test_try("Void constructor");
  1354 + try {
  1355 + GPhases phases;
  1356 + test_try_success();
  1357 + }
  1358 + catch (std::exception &e) {
  1359 + test_try_failure(e);
  1360 + }
  1361 +
  1362 + // Manipulate GPhases starting from an empty object
  1363 + GPhases phases;
  1364 + test_value(phases.size(), 0);
  1365 + test_assert(phases.is_empty(), "Test for empty GPhases object");
  1366 +
  1367 + // Add one phase interval
  1368 + phases.append(0.3, 0.5);
  1369 + test_value(phases.size(), 1);
  1370 + test_value(phases.pmin(0), 0.3);
  1371 + test_value(phases.pmax(0), 0.5);
  1372 + test_assert(!phases.is_empty(), "Test for non-empty GPhases object");
  1373 +
  1374 + // Remove phase interval
  1375 + phases.remove(0);
  1376 + test_value(phases.size(), 0);
  1377 + test_assert(phases.is_empty(), "Test for empty GPhases object");
  1378 +
  1379 + // Append two phase intervals
  1380 + phases.append(0.3, 0.5);
  1381 + phases.append(0.4, 0.7);
  1382 + test_value(phases.pmin(0), 0.3);
  1383 + test_value(phases.pmax(0), 0.5);
  1384 + test_value(phases.pmin(1), 0.4);
  1385 + test_value(phases.pmax(1), 0.7);
  1386 + test_value(phases.size(), 2);
  1387 + test_assert(!phases.is_empty(), "Test for non-empty GPhases object");
  1388 +
  1389 + // Copy phase intervals
  1390 + GPhases phases2(phases);
  1391 + test_value(phases2.pmin(0), 0.3);
  1392 + test_value(phases2.pmax(0), 0.5);
  1393 + test_value(phases2.pmin(1), 0.4);
  1394 + test_value(phases2.pmax(1), 0.7);
  1395 + test_value(phases2.size(), 2);
  1396 + test_assert(!phases2.is_empty(), "Test for non-empty GPhases object");
  1397 +
  1398 + // Clear object
  1399 + phases.clear();
  1400 + test_value(phases.size(), 0);
  1401 + test_assert(phases.is_empty(), "Test for empty GPhases object");
  1402 +
  1403 + // Extend photons
  1404 + phases.append(0.3, 0.5);
  1405 + phases.append(0.4, 0.7);
  1406 + phases.extend(phases);
  1407 + test_value(phases.size(), 4);
  1408 + test_assert(!phases.is_empty(), "Test for non-empty GPhases object");
  1409 +
  1410 + // Return
  1411 + return;
  1412 +}
  1413 +
  1414 +
  1415 +/***********************************************************************//**
1346 1416 * @brief Test GObservations::likelihood
1347 1417 ***************************************************************************/
1348 1418 void TestGObservation::test_observations_optimizer(void)
... ...
test/test_GObservation.hpp
  View file @9569954
... ... @@ -49,6 +49,7 @@ public:
49 49 virtual std::string classname(void) const { return "TestGObservation"; }
50 50 void test_ebounds(void);
51 51 void test_gti(void);
  52 + void test_phases(void);
52 53 void test_photons(void);
53 54 void test_time_reference(void);
54 55 void test_time(void);
... ...