diff --git a/ChangeLog b/ChangeLog index b674a32..91ffb63 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2017-04-11 +2017-04-14 * Version 1.3.0 released ========================= + Add GPhases class (#1996) Add GEventList::has_phase() and GEventList::has_detxy() methods (#1982) Add GEventList::append_column() method (#1991) Add scalar addition and subtraction operators to matrix classes (#1990) diff --git a/NEWS b/NEWS index 3947dac..f74e203 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.3.0.dev1 -11 April 2017 +14 April 2017 1. Introduction @@ -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 @@ -134,6 +135,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/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);