diff --git a/ChangeLog b/ChangeLog index 203e17a..c505a19 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.3.0 released ========================= + Add GPhases class (#1996) Add support for Monte Carlo identifier metadata (#1483) Add GEventList::has_mc_id() and GEventAtom::mc_id() methods (#1483) Add GEventList::has_phase() and GEventList::has_detxy() methods (#1982) diff --git a/NEWS b/NEWS index 9a32f9a..d3740d4 100644 --- a/NEWS +++ b/NEWS @@ -12,6 +12,7 @@ made since the last release of GammaLib. 2. Public interface modifications --------------------------------- The following classes have been added: +- GPhases - GWcsGLS - GWcsSFL - GModelSpectralSmoothBrokenPlaw @@ -139,6 +140,8 @@ GEbounds::set_log() method throws an exception if either the minimum or maximum energy are zero. Furthermore, the methods clear the energy boundaries if the requested number of intervals is non-zero (#1989). +The GPhases class has been added to handle phase selections of events (#1996). + 11. Optimization module ----------------------- diff --git a/include/GPhases.hpp b/include/GPhases.hpp new file mode 100644 index 0000000..ee8d66c --- /dev/null +++ b/include/GPhases.hpp @@ -0,0 +1,117 @@ +/*************************************************************************** + * GPhases.hpp - Phase intervals class * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GPhases.hpp + * @brief Phase intervals class interface definition + * @author Juergen Knoedlseder + */ + +#ifndef GPHASES_HPP +#define GPHASES_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GContainer.hpp" + +/* __ Forward declarations _______________________________________________ */ + + +/***********************************************************************//** + * @class GPhases + * + * @brief Phase Intervals class + ***************************************************************************/ +class GPhases : public GContainer { + +public: + // Constructors and destructors + GPhases(void); + GPhases(const GPhases& phases); + GPhases(const double& pmin, const double& pmax); + virtual ~GPhases(void); + + // Operators + GPhases& operator=(const GPhases& phases); + + // Methods + void clear(void); + GPhases* clone(void) const; + std::string classname(void) const; + int size(void) const; + bool is_empty(void) const; + void append(const double& pmin, const double& pmax); + void remove(const int& index); + void reserve(const int& num); + void extend(const GPhases& phases); + double pmin(const int& index) const; + double pmax(const int& index) const; + std::string print(const GChatter& chatter = NORMAL) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const GPhases& phases); + void free_members(void); + void insert_interval(const int& index, const double& pmin, + const double& pmax); + + // Protected data area + std::vector m_pmin; // Lower phase interval boundaries + std::vector m_pmax; // Upper phase interval boundaries +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GPhases"). + ***************************************************************************/ +inline +std::string GPhases::classname(void) const +{ + return ("GPhases"); +} + + +/***********************************************************************//** + * @brief Return number of phase intervals + * + * @return Number of phase intervals. + ***************************************************************************/ +inline +int GPhases::size(void) const +{ + return (int)m_pmin.size(); +} + + +/***********************************************************************//** + * @brief Signal if there are no phase intervals + * + * @return True if there are no phase intervals. + ***************************************************************************/ +inline +bool GPhases::is_empty(void) const +{ + return (m_pmin.size() == 0); +} + +#endif /* GPHASES_HPP */ diff --git a/include/GammaLib.hpp b/include/GammaLib.hpp index 5bf6bd7..bdfde50 100644 --- a/include/GammaLib.hpp +++ b/include/GammaLib.hpp @@ -183,6 +183,7 @@ #include "GGti.hpp" #include "GRoi.hpp" #include "GEbounds.hpp" +#include "GPhases.hpp" #include "GInstDir.hpp" #include "GResponse.hpp" #include "GPhotons.hpp" diff --git a/include/Makefile.am b/include/Makefile.am index 1c0f8df..f397498 100644 --- a/include/Makefile.am +++ b/include/Makefile.am @@ -143,6 +143,7 @@ pkginclude_HEADERS = GammaLib.hpp \ GGti.hpp \ GRoi.hpp \ GEbounds.hpp \ + GPhases.hpp \ GInstDir.hpp \ GResponse.hpp \ GPhotons.hpp \ diff --git a/inst/cta/include/GCTAEventList.hpp b/inst/cta/include/GCTAEventList.hpp index a0e4345..356d7d8 100644 --- a/inst/cta/include/GCTAEventList.hpp +++ b/inst/cta/include/GCTAEventList.hpp @@ -31,6 +31,7 @@ #include #include #include "GEventList.hpp" +#include "GPhases.hpp" #include "GFilename.hpp" #include "GCTAEventAtom.hpp" #include "GCTARoi.hpp" @@ -96,6 +97,8 @@ public: double irf_cache(const std::string& name, const int& index) const; void irf_cache(const std::string& name, const int& index, const double& irf) const; + const GPhases& phases(void) const; + void phases(const GPhases& phases); const std::string& gtiname(void) const; void has_phase(const bool& has_phase); void has_detxy(const bool& has_detxy); @@ -124,6 +127,7 @@ protected: // Event list meta data GCTARoi m_roi; //!< Region of interest + GPhases m_phases; //!< Phase intervals int m_num_events; //!< Number of events std::string m_gti_extname; //!< GTI extension name bool m_has_phase; //!< Signal presence of phase @@ -192,6 +196,31 @@ const GCTARoi& GCTAEventList::roi(void) const /***********************************************************************//** + * @brief Return phase intervals + * + * @return Phase intervals. + ***************************************************************************/ +inline +const GPhases& GCTAEventList::phases(void) const +{ + return (m_phases); +} + + +/***********************************************************************//** + * @brief Set phase intervals + * + * @param[in] phases Phase intervals. + ***************************************************************************/ +inline +void GCTAEventList::phases(const GPhases& phases) +{ + m_phases = phases; + return; +} + + +/***********************************************************************//** * @brief Reserves space for events * * @param[in] number Number of events. diff --git a/inst/cta/pyext/GCTAEventList.i b/inst/cta/pyext/GCTAEventList.i index cbf4876..60e5df0 100644 --- a/inst/cta/pyext/GCTAEventList.i +++ b/inst/cta/pyext/GCTAEventList.i @@ -62,17 +62,21 @@ public: void append_column(const GFitsTableCol& column); void reserve(const int& number); void remove(const int& index, const int& number = 1); - void write(GFits& fits, const std::string& evtname, - const std::string& gtiname) const; - void fetch(void); - void dispose(void); - double irf_cache(const std::string& name, const int& index) const; - void irf_cache(const std::string& name, const int& index, + void write(GFits& fits, + const std::string& evtname, + const std::string& gtiname) const; + void fetch(void) const; + void dispose(void) const; + double irf_cache(const std::string& name, + const int& index) const; + void irf_cache(const std::string& name, + const int& index, const double& irf) const; + const GPhases& phases(void) const; + void phases(const GPhases& phases); const std::string& gtiname(void) const; void has_phase(const bool& has_phase); void has_detxy(const bool& has_detxy); - void has_mc_id(const bool& has_mc_id); const bool& has_phase() const; const bool& has_detxy() const; const bool& has_mc_id() const; diff --git a/inst/cta/src/GCTAEventList.cpp b/inst/cta/src/GCTAEventList.cpp index 3fb6e11..2d2de0f 100644 --- a/inst/cta/src/GCTAEventList.cpp +++ b/inst/cta/src/GCTAEventList.cpp @@ -53,6 +53,8 @@ #define G_SET_MC_ID_NAMES "GCTAEventList::set_mc_id_names(std::vector&,"\ " std::vector&)" #define G_FETCH "GCTAEventList::fetch()" +#define G_WRITE_DS_KEYS "GCTAEventList::write_ds_keys(GFitsHDU&, "\ + "const std::string&) const" /* __ Macros _____________________________________________________________ */ @@ -439,6 +441,9 @@ void GCTAEventList::read(const GFits& fits) // Read energy boundaries from data sub-space keyword m_ebounds = gammalib::read_ds_ebounds(events); + // Read phase intervals from data sub-space keyword + m_phases = gammalib::read_ds_phase(events); + // If the file name contains an expression then load the events now. // We do this at the end to make sure that the GTI has been set before. if (m_filename.has_expression()) { @@ -957,6 +962,7 @@ void GCTAEventList::init_members(void) { // Initialise members m_roi.clear(); + m_phases.clear(); m_events.clear(); m_columns.clear(); m_filename.clear(); @@ -986,6 +992,7 @@ void GCTAEventList::copy_members(const GCTAEventList& list) { // Copy members m_roi = list.m_roi; + m_phases = list.m_phases; m_events = list.m_events; m_filename = list.m_filename; m_num_events = list.m_num_events; @@ -1337,6 +1344,9 @@ void GCTAEventList::write_events(GFitsBinTable& hdu) const * @param[in] hdu FITS HDU. * @param[in] gtiname Good Time Interval FITS extension. * + * @exception GException::vector_mismatch + * Lower and upper phase boundaries have different size. + * * Writes the data sub-space keywords for an event list into the FITS HDU. * The following keywords will be written: * @@ -1410,6 +1420,30 @@ void GCTAEventList::write_ds_keys(GFitsHDU& hdu, const std::string& gtiname) con } // endif: RoI was valid + // Check if there are phase interval cuts + if (m_phases.size() > 0) { + + // Set phase interval string + std::string dsval4 = ""; + for (int i = 0; i < m_phases.size(); ++i) { + if (i != 0) { + dsval4 += ","; + } + dsval4 += gammalib::str(m_phases.pmin(i)) + ":" + + gammalib::str(m_phases.pmax(i)); + } + + // Write DS keywords + hdu.card("DSTYP"+gammalib::str(ndskeys), "PHASE", + "Data sub-space type"); + hdu.card("DSUNI"+gammalib::str(ndskeys), "DIMENSIONLESS", + "Data sub-space unit"); + hdu.card("DSVAL"+gammalib::str(ndskeys), dsval4, + "Data sub-space value"); + ndskeys++; + + } // endif: there were phase interval cuts + // Set number of data selection keys hdu.card("NDSKEYS", ndskeys, "Number of data sub-space keys"); diff --git a/inst/cta/src/GCTASupport.cpp b/inst/cta/src/GCTASupport.cpp index 52b7dad..d13bb48 100644 --- a/inst/cta/src/GCTASupport.cpp +++ b/inst/cta/src/GCTASupport.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTASupport.cpp - CTA support functions * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2016 by Juergen Knoedlseder * + * copyright (C) 2011-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -36,13 +36,15 @@ #include "GMath.hpp" #include "GFitsHDU.hpp" #include "GEbounds.hpp" +#include "GPhases.hpp" #include "GCTARoi.hpp" #include "GCTAInstDir.hpp" /* __ Method name definitions ____________________________________________ */ #define G_READ_DS_ROI "gammalib::read_ds_roi(GFitsHDU&)" #define G_READ_DS_EBOUNDS "gammalib::read_ds_ebounds(GFitsHDU&)" -#define G_READ_DS_GTI "GGti gammalib::read_ds_gti_extname(GFitsHDU&)" +#define G_READ_DS_PHASE "gammalib::read_ds_phase(GFitsHDU&)" +#define G_READ_DS_GTI "gammalib::read_ds_gti_extname(GFitsHDU&)" /* __ Coding definitions _________________________________________________ */ @@ -259,6 +261,69 @@ GEbounds gammalib::read_ds_ebounds(const GFitsHDU& hdu) /***********************************************************************//** + * @brief Read phase boundary data sub-space keywords + * + * @param[in] hdu FITS HDU + * @return Phase intervals + * + * @exception GException::invalid_value + * Invalid phase data sub-space encountered + * + * Reads the phase boundary data sub-space keywords by searching for a + * DSTYPx keyword named "PHASE". The data sub-space information is expected + * to be in the format "0.1:0.3,0.5:0.7", where each subset of numbers + * represents the minimum and maximum phase. No detailed syntax + * checking is performed. + ***************************************************************************/ +GPhases gammalib::read_ds_phase(const GFitsHDU& hdu) +{ + // Initialise phase intervals + GPhases phases; + + // Get number of data sub-space keywords (default to 0 if keyword is + // not found) + int ndskeys = (hdu.has_card("NDSKEYS")) ? hdu.integer("NDSKEYS") : 0; + + // Loop over all data sub-space keys + for (int i = 1; i <= ndskeys; ++i) { + + // Set data sub-space key strings + std::string type_key = "DSTYP"+gammalib::str(i); + std::string unit_key = "DSUNI"+gammalib::str(i); + std::string value_key = "DSVAL"+gammalib::str(i); + + // Continue only if type_key is found and if this key is PHASE + if (hdu.has_card(type_key) && hdu.string(type_key) == "PHASE") { + + // Extract phase boundaries + std::string value = hdu.string(value_key); + std::vector intervals = gammalib::split(value, ","); + for (int j = 0; j < intervals.size(); ++j) { + std::vector values = gammalib::split(intervals[j], ":"); + if (values.size() == 2) { + double phmin = gammalib::todouble(values[0]); + double phmax = gammalib::todouble(values[1]); + phases.append(phmin, phmax); + } + else { + std::string msg = "Invalid phase value \""+value+ + "\" encountered in data selection key \""+ + value_key+"\""; + throw GException::invalid_value(G_READ_DS_PHASE, msg); + } + + } //endfor: looped over phase intervals + + } // endif: PHASE type_key found + + } // endfor: looped over data selection keys + + // Return + return phases; +} + + +/***********************************************************************//** * @brief Return Good Time Intervals extension name from data sub-space keywords * * @param[in] hdu FITS HDU diff --git a/inst/cta/src/GCTASupport.hpp b/inst/cta/src/GCTASupport.hpp index bdb723c..67e780f 100644 --- a/inst/cta/src/GCTASupport.hpp +++ b/inst/cta/src/GCTASupport.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTASupport.hpp - CTA support functions * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2016 by Juergen Knoedlseder * + * copyright (C) 2011-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -29,6 +29,7 @@ /* __ Includes ___________________________________________________________ */ #include +#include /* __ Namespaces _________________________________________________________ */ @@ -38,6 +39,7 @@ class GFitsHDU; class GCTARoi; class GEbounds; +class GPhases; /* __ Prototypes _________________________________________________________ */ namespace gammalib { @@ -46,6 +48,7 @@ namespace gammalib { const double& roi, const double& cosroi); GCTARoi read_ds_roi(const GFitsHDU& hdu); GEbounds read_ds_ebounds(const GFitsHDU& hdu); + GPhases read_ds_phase(const GFitsHDU& hdu); std::string read_ds_gti_extname(const GFitsHDU& hdu); } diff --git a/pyext/GGti.i b/pyext/GGti.i index 91ab6ac..400b5c8 100644 --- a/pyext/GGti.i +++ b/pyext/GGti.i @@ -35,7 +35,7 @@ * * @brief Good Time Interval class ***************************************************************************/ -class GGti : public GBase { +class GGti : public GContainer { public: // Constructors and destructors diff --git a/pyext/GPhases.i b/pyext/GPhases.i new file mode 100644 index 0000000..34b3521 --- /dev/null +++ b/pyext/GPhases.i @@ -0,0 +1,72 @@ +/*************************************************************************** + * GPhases.i - Phase intervals class * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GPhases.i + * @brief Phase intervals class interface definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GPhases.hpp" +#include "GTools.hpp" +%} + + +/***********************************************************************//** + * @class GPhases + * + * @brief Phase intervals class + ***************************************************************************/ +class GPhases : public GContainer { + +public: + // Constructors and destructors + GPhases(void); + GPhases(const GPhases& phases); + GPhases(const double& pmin, const double& pmax); + virtual ~GPhases(void); + + // Methods + void clear(void); + GPhases* clone(void) const; + std::string classname(void) const; + int size(void) const; + bool is_empty(void) const; + void append(const double& pmin, const double& pmax); + void remove(const int& index); + void reserve(const int& num); + void extend(const GPhases& phases); + double pmin(const int& index) const; + double pmax(const int& index) const; +}; + + +/***********************************************************************//** + * @brief GPhases class extension + ***************************************************************************/ +%extend GPhases { + int __len__() { + return (self->size()); + } + GPhases copy() { + return (*self); + } +}; diff --git a/pyext/gammalib/obs.i b/pyext/gammalib/obs.i index 9613a7b..6b47bc2 100644 --- a/pyext/gammalib/obs.i +++ b/pyext/gammalib/obs.i @@ -1,7 +1,7 @@ /*************************************************************************** * obs.i - Observation module * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2016 by Juergen Knoedlseder * + * copyright (C) 2011-2017 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -171,6 +171,7 @@ %include "GEventBin.i" %include "GInstDir.i" %include "GRoi.i" +%include "GPhases.i" %include "GResponse.i" %include "GPhoton.i" %include "GPhotons.i" diff --git a/src/obs/GGti.cpp b/src/obs/GGti.cpp index 22d6ebf..7521df3 100644 --- a/src/obs/GGti.cpp +++ b/src/obs/GGti.cpp @@ -100,7 +100,7 @@ GGti::GGti(const GFilename& filename) * * @param[in] gti Good Time Intervals. * - * Constructs Good Time Intervals by coping other Good Time Intervals. + * Constructs Good Time Intervals by copying other Good Time Intervals. ***************************************************************************/ GGti::GGti(const GGti& gti) { diff --git a/src/obs/GPhases.cpp b/src/obs/GPhases.cpp new file mode 100644 index 0000000..a52e64c --- /dev/null +++ b/src/obs/GPhases.cpp @@ -0,0 +1,483 @@ +/*************************************************************************** + * GPhases.cpp - Phase intervals class * + * ----------------------------------------------------------------------- * + * copyright (C) 2017 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GPhases.cpp + * @brief Phase intervals class implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include "GException.hpp" +#include "GTools.hpp" +#include "GPhases.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_REMOVE "GPhases::remove(int&)" +#define G_PMIN "GPhases::pmin(int&)" +#define G_PMAX "GPhases::pmax(int&)" +#define G_INSERT_INTERVAL "GPhases::insert_interval(int&, double&, double&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + * + * Constructs empty phase intervals. + ***************************************************************************/ +GPhases::GPhases(void) +{ + // Initialise class members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] phases Phase intervals. + * + * Constructs phase intervals by copying other phase intervals. + ***************************************************************************/ +GPhases::GPhases(const GPhases& phases) +{ + // Initialise class members + init_members(); + + // Copy members + copy_members(phases); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Single phase interval constructor + * + * @param[in] pmin Lower boundary of phase interval. + * @param[in] pmax Upper boundary of phase interval. + * + * Constructs Good Time Intervals from a single phase interval, given by + * [p@ pmin, @p pmax]. + ***************************************************************************/ +GPhases::GPhases(const double& pmin, const double& pmax) +{ + // Initialise members + init_members(); + + // Append phase interval + append(pmin, pmax); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GPhases::~GPhases(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] phases Phase intervals. + * @return Phase intervals. + ***************************************************************************/ +GPhases& GPhases::operator=(const GPhases& phases) +{ + // Execute only if object is not identical + if (this != &phases) { + + // Free members + free_members(); + + // Initialise private members + init_members(); + + // Copy members + copy_members(phases); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear phase intervals + ***************************************************************************/ +void GPhases::clear(void) +{ + // Free members + free_members(); + + // Initialise private members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone phase intervals + * + * @return Pointer to deep copy of phase intervals. + ***************************************************************************/ +GPhases* GPhases::clone(void) const +{ + return new GPhases(*this); +} + + +/***********************************************************************//** + * @brief Append phase interval + * + * @param[in] pmin Lower boundary of phase interval. + * @param[in] pmax Upper boundary of phase interval. + * + * Appends a phase interval at the end of the container. + ***************************************************************************/ +void GPhases::append(const double& pmin, const double& pmax) +{ + // Insert phase interval at the end of the list + insert_interval(size(), pmin, pmax); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Remove phase interval + * + * @param[in] index Phase interval index (0,...,size()-1). + * + * Removes phase interval at @p index from the container. All intervals + * after the specified @p index are moved forward by one position. + ***************************************************************************/ +void GPhases::remove(const int& index) +{ + #if defined(G_RANGE_CHECK) + // If index is outside boundary then throw an error + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_REMOVE, "Phase interval", index, + size()); + } + #endif + + // Remove lower and upper boundaries + m_pmin.erase(m_pmin.begin() + index); + m_pmax.erase(m_pmax.begin() + index); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Reserve space for phase intervals + * + * @param[in] num Number of elements. + ***************************************************************************/ +void GPhases::reserve(const int& num) +{ + // Reserves memory + m_pmin.reserve(num); + m_pmax.reserve(num); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Append phase intervals + * + * @param[in] phases Phase intervals. + * + * Append phase intervals to the container. + ***************************************************************************/ +void GPhases::extend(const GPhases& phases) +{ + // Continue only if phase intervals are not empty + if (!phases.is_empty()) { + + // Get size. Note that we extract the size first to avoid an + // endless loop that arises when a container is appended to + // itself. + int num = phases.size(); + + // Resize the boundary vectors + reserve(size() + num); + + // Loop over all elements and append them to container + for (int i = 0; i < num; ++i) { + m_pmin.push_back(phases.pmin(i)); + m_pmax.push_back(phases.pmax(i)); + } + + } // endif: Phase intervals were not empty + + // Return + return; +} + + +/***********************************************************************//** + * @brief Returns lower boundary for a given phase interval + * + * @param[in] index Phase interval index (0,...,size()-1). + * @return Lower boundary of phase interval. + * + * @exception GException::out_of_range + * Specified index is out of range. + ***************************************************************************/ +double GPhases::pmin(const int& index) const +{ + #if defined(G_RANGE_CHECK) + // If index is outside boundary then throw an error + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_PMIN, "Phase interval", index, + size()); + } + #endif + + // Return + return (m_pmin[index]); +} + + +/***********************************************************************//** + * @brief Returns upper boundary for a given phase interval + * + * @param[in] index Phase interval index (0,...,size()-1). + * @return Upper boundary of phase interval. + * + * @exception GException::out_of_range + * Specified index is out of range. + ***************************************************************************/ +double GPhases::pmax(const int& index) const +{ + #if defined(G_RANGE_CHECK) + // If index is outside boundary then throw an error + if (index < 0 || index >= size()) { + throw GException::out_of_range(G_PMAX, "Phase interval", index, + size()); + } + #endif + + // Return + return (m_pmax[index]); +} + + +/***********************************************************************//** + * @brief Print phase intervals + * + * @param[in] chatter Chattiness. + * @return String containing phase interval information. + ***************************************************************************/ +std::string GPhases::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GPhases ==="); + + // Append phase interval information + result.append("\n"+gammalib::parformat("Number of intervals")); + result.append(gammalib::str(size())); + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GPhases::init_members(void) +{ + // Initialise members + m_pmin.clear(); + m_pmax.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] phases Phase intervals. + ***************************************************************************/ +void GPhases::copy_members(const GPhases& phases) +{ + // Copy attributes + m_pmin = phases.m_pmin; + m_pmax = phases.m_pmax; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GPhases::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Insert phase interval + * + * @param[in] index Index after which interval is inserted. + * @param[in] pmin Lower boundary of interval. + * @param[in] pmax Upper boundary of interval. + * + * @exception GException::invalid_argument + * Invalid phase interval boundaries specified + * + * Inserts a phase interval before the interval with the specified @p index. + * If no interval with the specified index exists then append the interval + * at the end of the existing phase intervals. + ***************************************************************************/ +void GPhases::insert_interval(const int& index, + const double& pmin, + const double& pmax) +{ + // Throw an exception if phase minimum is outside the range [0,1] + if (pmin < 0.0 || pmin > 1.0) { + std::string msg = "Phase minimum "+gammalib::str(pmin)+ + " outside the valid range [0,1]. Please " + "specify phase interval boundaries comprised " + "within 0 and 1."; + throw GException::invalid_argument(G_INSERT_INTERVAL, msg); + } + + // Throw an exception if phase minimum is outside the range [0,1] + if (pmax < 0.0 || pmax > 1.0) { + std::string msg = "Phase maximum "+gammalib::str(pmax)+ + " outside the valid range [0,1]. Please " + "specify phase interval boundaries comprised " + "within 0 and 1."; + throw GException::invalid_argument(G_INSERT_INTERVAL, msg); + } + + // If phase minimum is smaller than phase maximum then append or insert + // a single interval ... + if (pmin < pmax) { + if (index >= size()) { + m_pmin.push_back(pmin); + m_pmax.push_back(pmax); + } + else { + m_pmin.insert(m_pmin.begin()+index, pmin); + m_pmax.insert(m_pmax.begin()+index, pmax); + } + } + + // ... otherwise if the phase minimum is larger than the maximum then + // consider this a wrap around interval and append or insert one interval + // from phase minimum to 1 and another interval from 0 to phase maximum + else if (pmin > pmax) { + if (pmin < 1.0) { + if (index >= size()) { + m_pmin.push_back(pmin); + m_pmax.push_back(1.0); + } + else { + m_pmin.insert(m_pmin.begin()+index, pmin); + m_pmax.insert(m_pmax.begin()+index, 1.0); + } + } + if (pmax > 0.0) { + if (index >= size()) { + m_pmin.push_back(0.0); + m_pmax.push_back(pmax); + } + else { + m_pmin.insert(m_pmin.begin()+index, 0.0); + m_pmax.insert(m_pmax.begin()+index, pmax); + } + } + } + + // Return + return; +} diff --git a/src/obs/Makefile.am b/src/obs/Makefile.am index dde7c05..0041332 100644 --- a/src/obs/Makefile.am +++ b/src/obs/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for obs module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016 Juergen Knoedlseder # +# Copyright (C) 2016-2017 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -44,6 +44,7 @@ sources = GEnergies.cpp \ GGti.cpp \ GRoi.cpp \ GEbounds.cpp \ + GPhases.cpp \ GResponse.cpp \ GInstDir.cpp \ GPhotons.cpp \ diff --git a/test/test_GObservation.cpp b/test/test_GObservation.cpp index 7324598..b68968d 100644 --- a/test/test_GObservation.cpp +++ b/test/test_GObservation.cpp @@ -69,6 +69,8 @@ void TestGObservation::set(void) "Test GEnergies class"); append(static_cast(&TestGObservation::test_ebounds), "Test GEbounds class"); + append(static_cast(&TestGObservation::test_phases), + "Test GPhases class"); append(static_cast(&TestGObservation::test_photons), "Test GPhotons class"); append(static_cast(&TestGObservation::test_observations_optimizer), @@ -1343,6 +1345,74 @@ void TestGObservation::test_photons(void) /***********************************************************************//** + * @brief Test GPhases + ***************************************************************************/ +void TestGObservation::test_phases(void) +{ + // Test void constructor + test_try("Void constructor"); + try { + GPhases phases; + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Manipulate GPhases starting from an empty object + GPhases phases; + test_value(phases.size(), 0); + test_assert(phases.is_empty(), "Test for empty GPhases object"); + + // Add one phase interval + phases.append(0.3, 0.5); + test_value(phases.size(), 1); + test_value(phases.pmin(0), 0.3); + test_value(phases.pmax(0), 0.5); + test_assert(!phases.is_empty(), "Test for non-empty GPhases object"); + + // Remove phase interval + phases.remove(0); + test_value(phases.size(), 0); + test_assert(phases.is_empty(), "Test for empty GPhases object"); + + // Append two phase intervals + phases.append(0.3, 0.5); + phases.append(0.4, 0.7); + test_value(phases.pmin(0), 0.3); + test_value(phases.pmax(0), 0.5); + test_value(phases.pmin(1), 0.4); + test_value(phases.pmax(1), 0.7); + test_value(phases.size(), 2); + test_assert(!phases.is_empty(), "Test for non-empty GPhases object"); + + // Copy phase intervals + GPhases phases2(phases); + test_value(phases2.pmin(0), 0.3); + test_value(phases2.pmax(0), 0.5); + test_value(phases2.pmin(1), 0.4); + test_value(phases2.pmax(1), 0.7); + test_value(phases2.size(), 2); + test_assert(!phases2.is_empty(), "Test for non-empty GPhases object"); + + // Clear object + phases.clear(); + test_value(phases.size(), 0); + test_assert(phases.is_empty(), "Test for empty GPhases object"); + + // Extend photons + phases.append(0.3, 0.5); + phases.append(0.4, 0.7); + phases.extend(phases); + test_value(phases.size(), 4); + test_assert(!phases.is_empty(), "Test for non-empty GPhases object"); + + // Return + return; +} + + +/***********************************************************************//** * @brief Test GObservations::likelihood ***************************************************************************/ void TestGObservation::test_observations_optimizer(void) diff --git a/test/test_GObservation.hpp b/test/test_GObservation.hpp index b6d96c9..9745920 100644 --- a/test/test_GObservation.hpp +++ b/test/test_GObservation.hpp @@ -49,6 +49,7 @@ public: virtual std::string classname(void) const { return "TestGObservation"; } void test_ebounds(void); void test_gti(void); + void test_phases(void); void test_photons(void); void test_time_reference(void); void test_time(void);