Eirik Hatlen / gammalib

Commit bed2fc7f812ebdf4ac41f604652b0b16ed52b62e

Authored by Jürgen Knödlseder
1 parent 6bf86700

Add CTA on-off observation handling classes. 

The GCTAOnOffObservation contains all information for a given run that is needed to perform an on/off analysis. The GCTAOnOffObservations container is a container class for GCTAOnOffObservation objects.
Inline Side-by-side
Showing 8 changed files with 1964 additions and 0 deletions
inst/cta/Makefile.am
  View file @bed2fc7
... ... @@ -11,6 +11,8 @@ INCLUDES = -I$(top_srcdir)/include \
11 11 sources = src/GCTASupport.cpp \
12 12 src/GCTAException.cpp \
13 13 src/GCTAObservation.cpp \
  14 + src/GCTAOnOffObservation.cpp \
  15 + src/GCTAOnOffObservations.cpp \
14 16 src/GCTAEventList.cpp \
15 17 src/GCTAEventAtom.cpp \
16 18 src/GCTAEventCube.cpp \
... ... @@ -41,6 +43,8 @@ sources = src/GCTASupport.cpp \
41 43 # Define headers to be installed
42 44 pkginclude_HEADERS = include/GCTAException.hpp \
43 45 include/GCTAObservation.hpp \
  46 + include/GCTAOnOffObservation.hpp \
  47 + include/GCTAOnOffObservations.hpp \
44 48 include/GCTAEventList.hpp \
45 49 include/GCTAEventAtom.hpp \
46 50 include/GCTAEventCube.hpp \
... ...
inst/cta/include/GCTAOnOffObservation.hpp 0 → 100644
  View file @bed2fc7
  1 +/***************************************************************************
  2 + * GCTAOnOffObservation.hpp - CTA on-off observation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Chia-Chun Lu & Christoph Deil *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservation.hpp
  23 + * @brief CTA on-off observation class definition
  24 + * @author Chia-Chun Lu & Christoph Deil
  25 + */
  26 +
  27 +#ifndef GCTAONOFFOBSERVATION_HPP
  28 +#define GCTAONOFFOBSERVATION_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GBase.hpp"
  33 +#include "GPha.hpp"
  34 +#include "GArf.hpp"
  35 +#include "GRmf.hpp"
  36 +#include "GCTAEventList.hpp"
  37 +#include "GCTAEventAtom.hpp"
  38 +#include "GCTAObservation.hpp"
  39 +#include "GSkyRegions.hpp"
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class GCTAOnOffObservation
  44 + *
  45 + * @brief CTA on-off observation class
  46 + ***************************************************************************/
  47 +class GCTAOnOffObservation : public GBase {
  48 +
  49 +public:
  50 + // Constructors and destructors
  51 + GCTAOnOffObservation(void);
  52 + GCTAOnOffObservation(const GEbounds& ereco, const GSkyRegions& on,
  53 + const GSkyRegions& off);
  54 + GCTAOnOffObservation(const GCTAOnOffObservation& obs);
  55 + virtual ~GCTAOnOffObservation(void);
  56 +
  57 + // Operators
  58 + GCTAOnOffObservation& operator=(const GCTAOnOffObservation& obs);
  59 +
  60 + // Methods
  61 + void clear(void);
  62 + GCTAOnOffObservation* clone(void) const;
  63 + void name(const std::string& name) { m_name = name; }
  64 + void instrument(const std::string& instrument) { m_instrument = instrument; }
  65 + void id(const std::string& id) { m_id = id; }
  66 + void off_regions(const GSkyRegions& regions) {m_off_regions = regions;}
  67 + void on_regions(const GSkyRegions& regions) {m_on_regions = regions;}
  68 + const std::string& name(void) const { return m_name; }
  69 + const std::string& instrument(void) const { return m_instrument; }
  70 + const std::string& id(void) const { return m_id; }
  71 + const GPha& on_spec(void) const { return m_on_spec; }
  72 + const GPha& off_spec(void) const { return m_off_spec; }
  73 + const GArf& arf(void) { return m_arf; };
  74 + const GRmf& rmf(void) {return m_rmf;}
  75 + void fill(const GCTAObservation& obs);
  76 + void compute_response(const GCTAObservation& obs,
  77 + const GEbounds& etrue);
  78 + void read(const GXmlElement& xml);
  79 + void write(GXmlElement& xml) const;
  80 + std::string print(const GChatter& chatter = NORMAL) const;
  81 +
  82 +protected:
  83 + // Protected methods
  84 + void init_members(void);
  85 + void copy_members(const GCTAOnOffObservation& obs);
  86 + void free_members(void);
  87 + void compute_arf(const GCTAObservation& obs);
  88 + void compute_rmf(const GCTAObservation& obs, const GEbounds& etrue);
  89 +
  90 + // Protected data members
  91 + std::string m_name; //!< Name
  92 + std::string m_instrument; //!< Instrument name
  93 + std::string m_id; //!< Observation identifier
  94 + GPha m_on_spec;
  95 + GPha m_off_spec;
  96 + GArf m_arf;
  97 + GRmf m_rmf;
  98 + GSkyRegions m_on_regions;
  99 + GSkyRegions m_off_regions;
  100 +};
  101 +
  102 +#endif /* GCTAONOFFOBSERVATION_HPP */
... ...
inst/cta/include/GCTAOnOffObservations.hpp 0 → 100644
  View file @bed2fc7
  1 +/***************************************************************************
  2 + * GCTAOnOffObservations.hpp - CTA on-off observation container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservations.hpp
  23 + * @brief CTA on-off observation container class definition
  24 + * @author Pierrick Martin
  25 + */
  26 +
  27 +#ifndef GCTAONOFFOBSERVATIONS_HPP
  28 +#define GCTAONOFFOBSERVATIONS_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include <vector>
  33 +#include "GContainer.hpp"
  34 +#include "GCTAOnOffObservation.hpp"
  35 +#include "GModels.hpp"
  36 +
  37 +
  38 +/***********************************************************************//**
  39 + * @class GCTAOnOffObservations
  40 + *
  41 + * @brief ON/OFF Observation container class
  42 + *
  43 + * This class in a container of GCTAOnOffObservation objects. Still some
  44 + * work to do and things to be clarified...
  45 + ***************************************************************************/
  46 +class GCTAOnOffObservations : public GContainer {
  47 +
  48 +public:
  49 + // Constructors and destructors
  50 + GCTAOnOffObservations(void);
  51 + GCTAOnOffObservations(const GCTAOnOffObservations& obs);
  52 + explicit GCTAOnOffObservations(const std::string& filename);
  53 + virtual ~GCTAOnOffObservations(void);
  54 +
  55 + // Operators
  56 + GCTAOnOffObservations& operator=(const GCTAOnOffObservations& obs);
  57 + GCTAOnOffObservation* operator[](const int& index);
  58 + const GCTAOnOffObservation* operator[](const int& index) const;
  59 +
  60 + // Methods
  61 + void clear(void);
  62 + GCTAOnOffObservations* clone(void) const;
  63 + GCTAOnOffObservation* at(const int& index);
  64 + const GCTAOnOffObservation* at(const int& index) const;
  65 + int size(void) const;
  66 + bool isempty(void) const;
  67 + GCTAOnOffObservation* set(const int& index, const GCTAOnOffObservation& obs);
  68 + GCTAOnOffObservation* append(const GCTAOnOffObservation& obs);
  69 + GCTAOnOffObservation* insert(const int& index, const GCTAOnOffObservation& obs);
  70 + void remove(const int& index);
  71 + void reserve(const int& num);
  72 + void extend(const GCTAOnOffObservations& obs);
  73 + bool contains(const std::string& instrument,
  74 + const std::string& id) const;
  75 + void load(const std::string& filename);
  76 + void save(const std::string& filename) const;
  77 + void read(const GXml& xml);
  78 + void write(GXml& xml) const;
  79 + void models(const GModels& models);
  80 + void models(const std::string& filename);
  81 + const GModels& models(void) const;
  82 + std::string print(const GChatter& chatter = NORMAL) const;
  83 +
  84 + /*
  85 + // To do
  86 + Implement fit functions for both cases summed over runs or not ?
  87 + Implement significance or have it at ctool level ?
  88 + Implement flux points computation and fit of a spectral function ?
  89 + // Do we want that ? If yes, need to retrieve inline functions from GObservations...
  90 + double npred(void) const;
  91 + // Left aside at the moment
  92 + void optimize(GOptimizer& opt);
  93 + */
  94 +
  95 +protected:
  96 + // Protected methods
  97 + void init_members(void);
  98 + void copy_members(const GCTAOnOffObservations& obs);
  99 + void free_members(void);
  100 + int get_index(const std::string& instrument,
  101 + const std::string& id) const;
  102 +
  103 + // Protected members
  104 + std::vector<GCTAOnOffObservation*> m_obs; //!< List of observations
  105 + GModels m_models; //!< List of models
  106 +};
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Return pointer to observation
  111 + *
  112 + * @param[in] index Observation index [0,...,size()-1].
  113 + * @return Observation.
  114 + *
  115 + * Returns a pointer to the observation with specified @p index.
  116 + ***************************************************************************/
  117 +inline
  118 +GCTAOnOffObservation* GCTAOnOffObservations::operator[](const int& index)
  119 +{
  120 + return m_obs[index];
  121 +}
  122 +
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Return pointer to observation (const version)
  126 + *
  127 + * @param[in] index Observation index [0,...,size()-1].
  128 + *
  129 + * Returns a const pointer to the observation with specified @p index.
  130 + ***************************************************************************/
  131 +inline
  132 +const GCTAOnOffObservation* GCTAOnOffObservations::operator[](const int& index) const
  133 +{
  134 + return m_obs[index];
  135 +}
  136 +
  137 +
  138 +/***********************************************************************//**
  139 + * @brief Return number of observations in container
  140 + *
  141 + * @return Number of observations in container.
  142 + *
  143 + * Returns the number of observations in the observation container.
  144 + ***************************************************************************/
  145 +inline
  146 +int GCTAOnOffObservations::size(void) const
  147 +{
  148 + return (m_obs.size());
  149 +}
  150 +
  151 +
  152 +/***********************************************************************//**
  153 + * @brief Signals if there are no observations in container
  154 + *
  155 + * @return True if container is empty, false otherwise.
  156 + *
  157 + * Signals if the observation container does not contain any observation.
  158 + ***************************************************************************/
  159 +inline
  160 +bool GCTAOnOffObservations::isempty(void) const
  161 +{
  162 + return (m_obs.empty());
  163 +}
  164 +
  165 +
  166 +/***********************************************************************//**
  167 + * @brief Reserves space for observations in container
  168 + *
  169 + * @param[in] num Number of observations.
  170 + *
  171 + * Reserves space for @p num observations in the container.
  172 + ***************************************************************************/
  173 +inline
  174 +void GCTAOnOffObservations::reserve(const int& num)
  175 +{
  176 + m_obs.reserve(num);
  177 + return;
  178 +}
  179 +
  180 +
  181 +/***********************************************************************//**
  182 + * @brief Set model container
  183 + *
  184 + * @param[in] models Model container.
  185 + *
  186 + * Sets the model container for the observations.
  187 + ***************************************************************************/
  188 +inline
  189 +void GCTAOnOffObservations::models(const GModels& models)
  190 +{
  191 + m_models = models;
  192 + return;
  193 +}
  194 +
  195 +
  196 +/***********************************************************************//**
  197 + * @brief Return model container
  198 + *
  199 + * @return Model container.
  200 + *
  201 + * Returns the model container of the observations.
  202 + ***************************************************************************/
  203 +inline
  204 +const GModels& GCTAOnOffObservations::models(void) const
  205 +{
  206 + return m_models;
  207 +}
  208 +
  209 +#endif /* GCTAONOFFOBSERVATIONS_HPP */
... ...
inst/cta/pyext/GCTAOnOffObservation.i 0 → 100644
  View file @bed2fc7
  1 +/***************************************************************************
  2 + * GCTAOnOffObservation.i - CTA on-off observation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservation.i
  23 + * @brief CTA on-off observation class definition
  24 + * @author Michael Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCTAOnOffObservation.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCTAOnOffObservation
  34 + *
  35 + * @brief CTA on-off observation class
  36 + ***************************************************************************/
  37 +class GCTAOnOffObservation : public GBase {
  38 +public:
  39 + // Constructors and destructors
  40 + GCTAOnOffObservation(void);
  41 + GCTAOnOffObservation(const GEbounds& ereco, const GSkyRegions& on,
  42 + const GSkyRegions& off);
  43 + GCTAOnOffObservation(const GCTAOnOffObservation& obs);
  44 + virtual ~GCTAOnOffObservation(void);
  45 +
  46 + // Methods
  47 + void clear(void);
  48 + GCTAOnOffObservation* clone(void) const;
  49 + void name(const std::string& name) { m_name = name; }
  50 + void instrument(const std::string& instrument) { m_instrument = instrument; }
  51 + void id(const std::string& id) { m_id = id; }
  52 + void off_regions(const GSkyRegions& regions) {m_off_regions = regions;}
  53 + void on_regions(const GSkyRegions& regions) {m_on_regions = regions;}
  54 + const std::string& name(void) const { return m_name; }
  55 + const std::string& instrument(void) const { return m_instrument; }
  56 + const std::string& id(void) const { return m_id; }
  57 + const GPha& on_spec(void) const { return m_on_spec; }
  58 + const GPha& off_spec(void) const { return m_off_spec; }
  59 + const GArf& arf(void) { return m_arf; };
  60 + const GRmf& rmf(void) {return m_rmf;}
  61 + void fill(const GCTAObservation& obs);
  62 + void compute_response(const GCTAObservation& obs,
  63 + const GEbounds& etrue);
  64 + void read(const GXmlElement& xml);
  65 + void write(GXmlElement& xml) const;
  66 +};
  67 +
  68 +
  69 +/***********************************************************************//**
  70 + * @brief GCTOnOffAObservation class extension
  71 + ***************************************************************************/
  72 +%extend GCTAOnOffObservation {
  73 + GCTAOnOffObservation copy() {
  74 + return (*self);
  75 + }
  76 +}
... ...
inst/cta/pyext/GCTAOnOffObservations.i 0 → 100644
  View file @bed2fc7
  1 +/***************************************************************************
  2 + * GCTAOnOffObservations.i - CTA on-off observation container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservations.i
  23 + * @brief CTA on-off observation container class definition
  24 + * @author Pierrick Martin
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCTAOnOffObservations.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCTAOnOffObservations
  34 + *
  35 + * @brief ON/OFF Observation container class
  36 + ***************************************************************************/
  37 +class GCTAOnOffObservations : public GContainer {
  38 +public:
  39 + // Constructors and destructors
  40 + GCTAOnOffObservations(void);
  41 + GCTAOnOffObservations(const GCTAOnOffObservations& obs);
  42 + explicit GCTAOnOffObservations(const std::string& filename);
  43 + virtual ~GCTAOnOffObservations(void);
  44 +
  45 + // Methods
  46 + void clear(void);
  47 + GCTAOnOffObservations* clone(void) const;
  48 + GCTAOnOffObservation* at(const int& index);
  49 + const GCTAOnOffObservation* at(const int& index) const;
  50 + int size(void) const;
  51 + bool isempty(void) const;
  52 + GCTAOnOffObservation* set(const int& index, const GCTAOnOffObservation& obs);
  53 + GCTAOnOffObservation* append(const GCTAOnOffObservation& obs);
  54 + GCTAOnOffObservation* insert(const int& index, const GCTAOnOffObservation& obs);
  55 + void remove(const int& index);
  56 + void reserve(const int& num);
  57 + void extend(const GCTAOnOffObservations& obs);
  58 + bool contains(const std::string& instrument,
  59 + const std::string& id) const;
  60 + void load(const std::string& filename);
  61 + void save(const std::string& filename) const;
  62 + void read(const GXml& xml);
  63 + void write(GXml& xml) const;
  64 + void models(const GModels& models);
  65 + void models(const std::string& filename);
  66 + const GModels& models(void) const;
  67 +};
  68 +
  69 +
  70 +/***********************************************************************//**
  71 + * @brief GCTAOnOffObservations class extension
  72 + ***************************************************************************/
  73 +%extend GCTAOnOffObservations {
  74 + GCTAOnOffObservation* __getitem__(const int& index) {
  75 + if (index >= 0 && index < self->size()) {
  76 + return (*self)[index];
  77 + }
  78 + else {
  79 + throw GException::out_of_range("__getitem__(int)", index, self->size());
  80 + }
  81 + }
  82 + void __setitem__(const int& index, const GCTAOnOffObservation& val) {
  83 + if (index>=0 && index < self->size()) {
  84 + self->set(index, val);
  85 + return;
  86 + }
  87 + else {
  88 + throw GException::out_of_range("__setitem__(int)", index, self->size());
  89 + }
  90 + }
  91 + int __len__() {
  92 + return (self->size());
  93 + }
  94 + GCTAOnOffObservations copy() {
  95 + return (*self);
  96 + }
  97 +};
... ...
inst/cta/pyext/cta.i
  View file @bed2fc7
... ... @@ -63,6 +63,8 @@
63 63  
64 64 /* __ CTA ________________________________________________________________ */
65 65 %include "GCTAObservation.i"
  66 +%include "GCTAOnOffObservation.i"
  67 +%include "GCTAOnOffObservations.i"
66 68 %include "GCTAEventCube.i"
67 69 %include "GCTAEventList.i"
68 70 %include "GCTAEventBin.i"
... ...
inst/cta/src/GCTAOnOffObservation.cpp 0 → 100644
  View file @bed2fc7
  1 +/***************************************************************************
  2 + * GCTAOnOffObservation.cpp - CTA on-off observation class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Chia-Chun Lu & Christoph Deil *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservation.cpp
  23 + * @brief CTA on-off observation class implementation
  24 + * @author Chia-Chun Lu & Christoph Deil
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GCTAOnOffObservation.hpp"
  32 +#include "GTools.hpp"
  33 +
  34 +/* __ Method name definitions ____________________________________________ */
  35 +#define G_WRITE "GCTAOnOffObservation::write(GXmlElement&)"
  36 +#define G_READ "GCTAOnOffObservation::read(GXmlElement&)"
  37 +#define G_FILL "GCTAOnOffObservation::fill(GCTAObservation&)"
  38 +
  39 +/* __ Macros _____________________________________________________________ */
  40 +
  41 +/* __ Coding definitions _________________________________________________ */
  42 +
  43 +/* __ Debug definitions __________________________________________________ */
  44 +
  45 +
  46 +/*==========================================================================
  47 + = =
  48 + = Constructors/destructors =
  49 + = =
  50 + ==========================================================================*/
  51 +
  52 +/***********************************************************************//**
  53 + * @brief Void constructor
  54 + ***************************************************************************/
  55 +GCTAOnOffObservation::GCTAOnOffObservation(void)
  56 +{
  57 + // Initialise private members
  58 + init_members();
  59 +
  60 + // Return
  61 + return;
  62 +}
  63 +
  64 +
  65 +/***********************************************************************//**
  66 + * @brief Copy constructor
  67 + *
  68 + * @param[in] obs CTA on-off observation.
  69 + ***************************************************************************/
  70 +GCTAOnOffObservation::GCTAOnOffObservation(const GCTAOnOffObservation& obs)
  71 +{
  72 + // Initialise private
  73 + init_members();
  74 +
  75 + // Copy members
  76 + copy_members(obs);
  77 +
  78 + // Return
  79 + return;
  80 +}
  81 +
  82 +
  83 +/***********************************************************************//**
  84 + * @brief Parameter constructor
  85 + *
  86 + * @param[in] ereco Reconstructed energy bins.
  87 + * @param[in] on On regions.
  88 + * @param[in] off Off regions.
  89 + ***************************************************************************/
  90 +GCTAOnOffObservation::GCTAOnOffObservation(const GEbounds& ereco,
  91 + const GSkyRegions& on,
  92 + const GSkyRegions& off)
  93 +{
  94 + // Initialise private
  95 + init_members();
  96 +
  97 + // Initialise spectra
  98 + m_on_spec = GPha(ereco);
  99 + m_off_spec = GPha(ereco);
  100 +
  101 + // Store regions
  102 + m_on_regions = on;
  103 + m_off_regions = off;
  104 +
  105 + // Return
  106 + return;
  107 +}
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Destructor
  112 + ***************************************************************************/
  113 +GCTAOnOffObservation::~GCTAOnOffObservation(void)
  114 +{
  115 + // Free members
  116 + free_members();
  117 +
  118 + // Return
  119 + return;
  120 +}
  121 +
  122 +
  123 +/*==========================================================================
  124 + = =
  125 + = Operators =
  126 + = =
  127 + ==========================================================================*/
  128 +
  129 +/***********************************************************************//**
  130 + * @brief Assignment operator
  131 + *
  132 + * @param[in] obs CTA on-off observation.
  133 + ***************************************************************************/
  134 +GCTAOnOffObservation& GCTAOnOffObservation::operator=(const GCTAOnOffObservation& obs)
  135 +{
  136 + // Execute only if object is not identical
  137 + if (this != &obs) {
  138 +
  139 + // Free members
  140 + free_members();
  141 +
  142 + // Initialise private members for clean destruction
  143 + init_members();
  144 +
  145 + // Copy members
  146 + copy_members(obs);
  147 +
  148 + } // endif: object was not identical
  149 +
  150 + // Return
  151 + return *this;
  152 +}
  153 +
  154 +
  155 +/*==========================================================================
  156 + = =
  157 + = Public methods =
  158 + = =
  159 + ==========================================================================*/
  160 +
  161 +/***********************************************************************//**
  162 + * @brief Clear instance
  163 + *
  164 + * This method properly resets the object to an initial state.
  165 + ***************************************************************************/
  166 +void GCTAOnOffObservation::clear(void)
  167 +{
  168 + // Free class members
  169 + free_members();
  170 +
  171 + // Initialise members
  172 + init_members();
  173 +
  174 + // Return
  175 + return;
  176 +}
  177 +
  178 +
  179 +/***********************************************************************//**
  180 + * @brief Clone instance
  181 + **************************************************************************/
  182 +GCTAOnOffObservation* GCTAOnOffObservation::clone(void) const
  183 +{
  184 + return new GCTAOnOffObservation(*this);
  185 +}
  186 +
  187 +
  188 +/***********************************************************************//**
  189 + * @brief read observation from an xml element
  190 + *
  191 + * @param[in] xml GXmlElement.
  192 + ***************************************************************************/
  193 +void GCTAOnOffObservation::read(const GXmlElement& xml)
  194 +{
  195 + // clean object
  196 + clear();
  197 +
  198 + // Extract instrument name
  199 + m_instrument = xml.attribute("instrument");
  200 +
  201 + // Determine number of parameter nodes in XML element
  202 + int npars = xml.elements("parameter");
  203 +
  204 + // Verify that XML element has exactly 6 parameters
  205 + if (xml.elements() != 6 || npars != 6) {
  206 + throw GException::xml_invalid_parnum(G_READ, xml,
  207 + "CTA observation requires exactly 6 parameters.");
  208 + }
  209 +
  210 + // Extract parameters
  211 + int npar[] = {0, 0, 0, 0, 0, 0};
  212 + for (int i = 0; i < npars; ++i) {
  213 +
  214 + // Get parameter element
  215 + const GXmlElement* par = xml.element("parameter", i);
  216 +
  217 + // Handle Pha_on
  218 + if (par->attribute("name") == "Pha_on") {
  219 +
  220 + // Read pha_on file name
  221 + std::string filename = par->attribute("file");
  222 +
  223 + // Load pha_on
  224 + m_on_spec.load(filename);
  225 +
  226 + // Increment parameter counter
  227 + npar[0]++;
  228 + }
  229 +
  230 + // Handle Pha_off
  231 + else if (par->attribute("name") == "Pha_off") {
  232 +
  233 + // Read countsmap file name
  234 + std::string filename = par->attribute("file");
  235 +
  236 + // Load pha_off
  237 + m_off_spec.load(filename);
  238 +
  239 + // Increment parameter counter
  240 + npar[1]++;
  241 + }
  242 +
  243 + // Handle on regions
  244 + else if(par->attribute("name") == "Regions_on") {
  245 +
  246 + // Get filename
  247 + std::string filename = par->attribute("file");
  248 +
  249 + // load on regions
  250 + m_on_regions.load(filename);
  251 +
  252 + // Increase number of parameters
  253 + npar[2]++;
  254 + }
  255 +
  256 + // Handle off regions
  257 + else if(par->attribute("name") == "Regions_off") {
  258 +
  259 + // Get filename
  260 + std::string filename = par->attribute("file");
  261 +
  262 + // load off regions
  263 + m_on_regions.load(filename);
  264 +
  265 + // Increase number of parameters
  266 + npar[3]++;
  267 + }
  268 +
  269 + // Handle Arf
  270 + else if(par->attribute("name") == "Arf") {
  271 +
  272 + // Get filename
  273 + std::string filename = par->attribute("file");
  274 +
  275 + // load arf
  276 + m_arf.load(filename);
  277 +
  278 + // Increase number of parameters
  279 + npar[4]++;
  280 + }
  281 +
  282 + // Handle Rmf
  283 + else if(par->attribute("name") == "Rmf") {
  284 +
  285 + // Get filename
  286 + std::string filename = par->attribute("file");
  287 +
  288 + // load Rmf
  289 + m_rmf.load(filename);
  290 +
  291 + // Increase number of parameters
  292 + npar[5]++;
  293 + }
  294 +
  295 + } // endfor: looped over all parameters
  296 +
  297 + // Verify that all parameters were found
  298 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1 || npar[4] != 1 || npar[5] != 1) {
  299 + throw GException::xml_invalid_parnames(G_READ, xml,
  300 + "Require \"Pha_on\" or \"Pha_off\" and \"Regions_on\""
  301 + ", \"Regions_off\",\"Rmf\" and \"Arf\" parameters.");
  302 + }
  303 +
  304 + // Return
  305 + return;
  306 +}
  307 +
  308 +
  309 +/***********************************************************************//**
  310 + * @brief write observation to an xml element
  311 + *
  312 + * @param[in] xml GXmlElement.
  313 + ***************************************************************************/
  314 +void GCTAOnOffObservation::write(GXmlElement& xml) const
  315 +{
  316 + // If XML element has 0 nodes then append 6 parameter nodes
  317 + if (xml.elements() == 0) {
  318 + xml.append(GXmlElement("parameter name=\"Pha_on\""));
  319 + xml.append(GXmlElement("parameter name=\"Pha_off\""));
  320 + xml.append(GXmlElement("parameter name=\"Regions_on\""));
  321 + xml.append(GXmlElement("parameter name=\"Regions_off\""));
  322 + xml.append(GXmlElement("parameter name=\"Arf\""));
  323 + xml.append(GXmlElement("parameter name=\"Rmf\""));
  324 + }
  325 +
  326 + // Verify that XML element has exactly 4 parameters
  327 + if (xml.elements() != 6 || xml.elements("parameter") != 6) {
  328 + throw GException::xml_invalid_parnum(G_WRITE, xml,
  329 + "CTAOnOffObservation requires exactly 6 parameters.");
  330 + }
  331 +
  332 + // Set or update parameter attributes
  333 + int npar[] = {0, 0, 0, 0, 0, 0};
  334 + for (int i = 0; i < 6; ++i) {
  335 +
  336 + // Get parameter element
  337 + GXmlElement* par = xml.element("parameter", i);
  338 +
  339 + // Handle on counts
  340 + if (par->attribute("name") == "Pha_on") {
  341 + par->attribute("file", m_on_spec.filename());
  342 + npar[0]++;
  343 + }
  344 +
  345 + // Handle off counts
  346 + else if (par->attribute("name") == "Pha_off") {
  347 + par->attribute("file", m_off_spec.filename());
  348 + npar[0]++;
  349 + }
  350 +
  351 + // Handle on regions
  352 + else if (par->attribute("name") == "Regions_on") {
  353 + par->attribute("file", m_on_regions.filename());
  354 + npar[0]++;
  355 + }
  356 +
  357 + // Handle off regions
  358 + else if (par->attribute("name") == "Regions_off") {
  359 + par->attribute("file", m_off_regions.filename());
  360 + npar[0]++;
  361 + }
  362 +
  363 + // Handle effective area
  364 + else if (par->attribute("name") == "Arf") {
  365 + par->attribute("file", m_arf.filename());
  366 + npar[0]++;
  367 + }
  368 +
  369 + // Handle energy resolution
  370 + else if (par->attribute("name") == "Rmf") {
  371 + par->attribute("file", m_rmf.filename());
  372 + npar[0]++;
  373 + }
  374 + }
  375 + // Verify that all required parameters are present
  376 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1 || npar[4] != 1 || npar[5] !=1) {
  377 + throw GException::xml_invalid_parnames(G_WRITE, xml,
  378 + "Require \"Pha_on\" or \"Pha_off\" and \"Regions_on\", \"Regions_off\""
  379 + " \"Arf\", \"Rmf\" parameters.");
  380 + }
  381 +
  382 + // Return
  383 + return;
  384 +}
  385 +
  386 +
  387 +/***********************************************************************//**
  388 + * @brief Fill events in the on and off spectra
  389 + *
  390 + * @param[in] obs CTA observation.
  391 + *
  392 + * @exception GException::invalid_value
  393 + * No CTA event list found in CTA observation.
  394 + ***************************************************************************/
  395 +void GCTAOnOffObservation::fill(const GCTAObservation& obs)
  396 +{
  397 + // Get CTA event list pointer
  398 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  399 + if (events == NULL) {
  400 + throw GException::invalid_value(G_FILL,
  401 + "GCTAObservation must provide an event list");
  402 + }
  403 +
  404 + // Loop over all events
  405 + for (int i = 0; i < events->size(); ++i) {
  406 +
  407 + // Get measured event direction
  408 + const GCTAEventAtom* atom = (*events)[i];
  409 + GSkyDir dir = atom->dir().dir();
  410 +
  411 + // Fill in spectrum according to region containment
  412 + if (m_on_regions.contains(dir)) {
  413 + m_on_spec.fill(atom->energy());
  414 + }
  415 + else if (m_off_regions.contains(dir)) {
  416 + m_off_spec.fill(atom->energy());
  417 + }
  418 +
  419 + } // endfor: looped over all events
  420 +
  421 + // Return
  422 + return;
  423 +}
  424 +
  425 +
  426 +/***********************************************************************//**
  427 + * @brief Compute response of an on/off observation
  428 + *
  429 + * @param[in] obs GCTAObservation.
  430 + * @param[in] etrue True energy boundaries.
  431 + ***************************************************************************/
  432 +void GCTAOnOffObservation::compute_response(const GCTAObservation& obs,
  433 + const GEbounds& etrue)
  434 +{
  435 + // Compute response components
  436 + compute_arf(obs);
  437 + compute_rmf(obs, etrue);
  438 +
  439 + // Return
  440 + return;
  441 +}
  442 +
  443 +/***********************************************************************//**
  444 + * @brief Compute arf of an on/off observation
  445 + *
  446 + * @param[in] obs GCTAObservation.
  447 + *
  448 + * @todo Implement GCTAResponse::npred usage.
  449 + ***************************************************************************/
  450 +void GCTAOnOffObservation::compute_arf(const GCTAObservation& obs)
  451 +{
  452 + // Set constant response parameters
  453 + const double theta = 0.0;
  454 + const double phi = 0.0;
  455 + const double zenith = 0.0;
  456 + const double azimuth = 0.0;
  457 +
  458 + // Get energy boundaries from on spectrum
  459 + GEbounds ereco = m_on_spec.ebounds();
  460 +
  461 + // Continue only if there are spectral bins
  462 + int nreco = ereco.size();
  463 + if (nreco > 0) {
  464 +
  465 + // Get CTA response pointer
  466 + GCTAResponse* response = obs.response();
  467 +
  468 + // Initialize ARF
  469 + m_arf = GArf(ereco);
  470 +
  471 + // Loop over reconstructed energies
  472 + for (int i = 0; i < nreco; ++i) {
  473 +
  474 + // Get mean energy of bin
  475 + double logenergy = ereco.elogmean(i).log10TeV();
  476 +
  477 + // Set specresp value
  478 + m_arf[i] = response->aeff(theta,
  479 + phi,
  480 + zenith,
  481 + azimuth,
  482 + logenergy);
  483 +
  484 + } // endfor: looped over reconstructed energies
  485 +
  486 + } // endif: there were energy bins
  487 +
  488 + // Return
  489 + return;
  490 +}
  491 +
  492 +
  493 +/***********************************************************************//**
  494 + * @brief compute rmf of an on/off observation
  495 + *
  496 + * @param[in] obs GCTAObservation.
  497 + * @param[in] etrue True energy boundaries.
  498 + ***************************************************************************/
  499 +void GCTAOnOffObservation::compute_rmf(const GCTAObservation& obs,
  500 + const GEbounds& etrue)
  501 +{
  502 + // Set constant response parameters
  503 + const double theta = 0.0;
  504 + const double phi = 0.0;
  505 + const double zenith = 0.0;
  506 + const double azimuth = 0.0;
  507 +
  508 + // Get reconstructed energy boundaries from on spectrum
  509 + GEbounds ereco = m_on_spec.ebounds();
  510 +
  511 + // Continue only if there are spectral bins
  512 + int ntrue = etrue.size();
  513 + int nreco = ereco.size();
  514 + if (ntrue > 0 && nreco > 0) {
  515 +
  516 + // Get CTA response pointer
  517 + GCTAResponse* response = obs.response();
  518 +
  519 + // Initialize RMF
  520 + m_rmf = GRmf(etrue, ereco);
  521 +
  522 + // Loop over reconstructed energy
  523 + for (int ireco = 0; ireco < nreco; ++ireco) {
  524 +
  525 + // Compute reconstructed energy
  526 + double eng_reco = ereco.elogmean(ireco).log10TeV();
  527 +
  528 + // Loop over true energy
  529 + for (int itrue = 0; itrue < ntrue; ++itrue) {
  530 +
  531 + // Compute true energy
  532 + double eng_true = etrue.elogmean(itrue).log10TeV();
  533 +
  534 + // Set RMF value
  535 + m_rmf(itrue, ireco) = response->edisp(eng_reco,
  536 + theta,
  537 + phi,
  538 + zenith,
  539 + azimuth,
  540 + eng_true);
  541 + } // endfor: looped over true energy
  542 + } // endfor: looped over reconstructed energy
  543 + } // endif: there were energy bins
  544 +
  545 + // Return
  546 + return;
  547 +}
  548 +
  549 +
  550 +/***********************************************************************//**
  551 + * @brief Print class information
  552 + ***************************************************************************/
  553 +std::string GCTAOnOffObservation::print(const GChatter& chatter) const
  554 +{
  555 + // Initialise result string
  556 + std::string result;
  557 +
  558 + // Continue only if chatter is not silent
  559 + if (chatter != SILENT) {
  560 +
  561 + // Append header
  562 + result.append("=== GCTAOnOffObservation ===");
  563 +
  564 + // Append parameters
  565 + result.append("\n"+gammalib::parformat("Name")+m_name);
  566 + result.append("\n"+gammalib::parformat("Identifier")+m_id);
  567 +
  568 + // Append spectra, ARF and RMF
  569 + result.append("\n"+m_on_spec.print());
  570 + result.append("\n"+m_off_spec.print());
  571 + result.append("\n"+m_arf.print());
  572 + result.append("\n"+m_rmf.print());
  573 +
  574 + // Append regions
  575 + result.append("\n"+m_on_regions.print());
  576 + result.append("\n"+m_off_regions.print());
  577 + }
  578 +
  579 + // Return result
  580 + return result;
  581 +}
  582 +
  583 +
  584 +/*==========================================================================
  585 + = =
  586 + = Private methods =
  587 + = =
  588 + ==========================================================================*/
  589 +
  590 +/***********************************************************************//**
  591 + * @brief Initialise class members
  592 + ***************************************************************************/
  593 +void GCTAOnOffObservation::init_members(void)
  594 +{
  595 + // Initialise members
  596 + m_name.clear();
  597 + m_id.clear();
  598 + m_instrument = "CTA";
  599 + m_on_spec.clear();
  600 + m_off_spec.clear();
  601 + m_arf.clear();
  602 + m_rmf.clear();
  603 + m_on_regions.clear();
  604 + m_off_regions.clear();
  605 +
  606 + // Return
  607 + return;
  608 +}
  609 +
  610 +
  611 +/***********************************************************************//**
  612 + * @brief Copy class members
  613 + *
  614 + * @param[in] obs CTA on-off observation.
  615 + ***************************************************************************/
  616 +void GCTAOnOffObservation::copy_members(const GCTAOnOffObservation& obs)
  617 +{
  618 + // Copy attributes
  619 + m_name = obs.m_name;
  620 + m_id = obs.m_id;
  621 + m_instrument = obs.m_instrument;
  622 + m_on_spec = obs.m_on_spec;
  623 + m_off_spec = obs.m_off_spec;
  624 + m_arf = obs.m_arf;
  625 + m_rmf = obs.m_rmf;
  626 + m_on_regions = obs.m_on_regions;
  627 + m_off_regions = obs.m_off_regions;
  628 +
  629 + // Return
  630 + return;
  631 +}
  632 +
  633 +
  634 +
  635 +/***********************************************************************//**
  636 + * @brief Delete class members
  637 + ***************************************************************************/
  638 +void GCTAOnOffObservation::free_members(void)
  639 +{
  640 + // Return
  641 + return;
  642 +}
... ...
inst/cta/src/GCTAOnOffObservations.cpp 0 → 100644
  View file @bed2fc7
  1 +/***************************************************************************
  2 + * GCTAOnOffObservations.cpp - ON/OFF Observation container class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Pierrick Martin *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAOnOffObservations.cpp
  23 + * @brief ON/OFF Observations container class implementation
  24 + * @author Pierrick Martin
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GException.hpp"
  33 +#include "GCTAOnOffObservations.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_AT "GCTAOnOffObservations::at(int&)"
  37 +#define G_SET "GCTAOnOffObservations::set(int&, GCTAOnOffObservation&)"
  38 +#define G_APPEND "GCTAOnOffObservations::append(GCTAOnOffObservation&)"
  39 +#define G_INSERT "GCTAOnOffObservations::insert(int&, GCTAOnOffObservation&)"
  40 +#define G_REMOVE "GCTAOnOffObservations::remove(int&)"
  41 +#define G_EXTEND "GCTAOnOffObservations::extend(GCTAOnOffObservations&)"
  42 +#define G_READ "GCTAOnOffObservations::read(GXml&)"
  43 +
  44 +/* __ Macros _____________________________________________________________ */
  45 +
  46 +/* __ Coding definitions _________________________________________________ */
  47 +
  48 +/* __ Debug definitions __________________________________________________ */
  49 +
  50 +
  51 +/*==========================================================================
  52 + = =
  53 + = Constructors/destructors =
  54 + = =
  55 + ==========================================================================*/
  56 +
  57 +/***********************************************************************//**
  58 + * @brief Void constructor
  59 + ***************************************************************************/
  60 +GCTAOnOffObservations::GCTAOnOffObservations(void)
  61 +{
  62 + // Initialise members
  63 + init_members();
  64 +
  65 + // Return
  66 + return;
  67 +}
  68 +
  69 +
  70 +/***********************************************************************//**
  71 + * @brief Copy constructor
  72 + *
  73 + * @param obs Observation container.
  74 + ***************************************************************************/
  75 +GCTAOnOffObservations::GCTAOnOffObservations(const GCTAOnOffObservations& obs)
  76 +{
  77 + // Initialise members
  78 + init_members();
  79 +
  80 + // Copy members
  81 + copy_members(obs);
  82 +
  83 + // Return
  84 + return;
  85 +}
  86 +
  87 +
  88 +/***********************************************************************//**
  89 + * @brief Load constructor
  90 + *
  91 + * @param[in] filename XML filename.
  92 + *
  93 + * Construct the observation container by loading all relevant information
  94 + * from a XML file. Please refer to the read() method for more information
  95 + * about the structure of the XML file.
  96 + ***************************************************************************/
  97 +GCTAOnOffObservations::GCTAOnOffObservations(const std::string& filename)
  98 +{
  99 + // Initialise members
  100 + init_members();
  101 +
  102 + // Load XML file
  103 + load(filename);
  104 +
  105 + // Return
  106 + return;
  107 +}
  108 +
  109 +
  110 +/***********************************************************************//**
  111 + * @brief Destructor
  112 + ***************************************************************************/
  113 +GCTAOnOffObservations::~GCTAOnOffObservations(void)
  114 +{
  115 + // Free members
  116 + free_members();
  117 +
  118 + // Return
  119 + return;
  120 +}
  121 +
  122 +
  123 +/*==========================================================================
  124 + = =
  125 + = Operators =
  126 + = =
  127 + ==========================================================================*/
  128 +
  129 +/***********************************************************************//**
  130 + * @brief Assignment operator
  131 + *
  132 + * @param[in] obs Observation container.
  133 + * @return Observation container.
  134 + ***************************************************************************/
  135 +GCTAOnOffObservations& GCTAOnOffObservations::operator=(const GCTAOnOffObservations& obs)
  136 +{
  137 + // Execute only if object is not identical
  138 + if (this != &obs) {
  139 +
  140 + // Free members
  141 + free_members();
  142 +
  143 + // Initialise members
  144 + init_members();
  145 +
  146 + // Copy members
  147 + copy_members(obs);
  148 +
  149 + } // endif: object was not identical
  150 +
  151 + // Return this object
  152 + return *this;
  153 +}
  154 +
  155 +
  156 +/*==========================================================================
  157 + = =
  158 + = Public methods =
  159 + = =
  160 + ==========================================================================*/
  161 +
  162 +/***********************************************************************//**
  163 + * @brief Clear observations
  164 + ***************************************************************************/
  165 +void GCTAOnOffObservations::clear(void)
  166 +{
  167 + // Free members
  168 + free_members();
  169 +
  170 + // Initialise members
  171 + init_members();
  172 +
  173 + // Return
  174 + return;
  175 +}
  176 +
  177 +
  178 +/***********************************************************************//**
  179 + * @brief Clone observations
  180 + *
  181 + * @return Pointer to deep copy of observation container.
  182 + ***************************************************************************/
  183 +GCTAOnOffObservations* GCTAOnOffObservations::clone(void) const
  184 +{
  185 + // Clone observations
  186 + return new GCTAOnOffObservations(*this);
  187 +}
  188 +
  189 +
  190 +/***********************************************************************//**
  191 + * @brief Return pointer to observation
  192 + *
  193 + * @param[in] index Observation index [0,...,size()-1].
  194 + * @return Observation.
  195 + *
  196 + * @exception GException::out_of_range
  197 + * Operation index is out of range.
  198 + *
  199 + * Returns a pointer to the observation with specified @p index.
  200 + ***************************************************************************/
  201 +GCTAOnOffObservation* GCTAOnOffObservations::at(const int& index)
  202 +{
  203 + // Raise exception if index is out of range
  204 + if (index < 0 || index >= size()) {
  205 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  206 + }
  207 +
  208 + // Return pointer
  209 + return m_obs[index];
  210 +}
  211 +
  212 +
  213 +/***********************************************************************//**
  214 + * @brief Return pointer to observation (const version)
  215 + *
  216 + * @param[in] index Observation index [0,...,size()-1].
  217 + *
  218 + * @exception GException::out_of_range
  219 + * Operation index is out of range.
  220 + *
  221 + * Returns a const pointer to the observation with specified @p index.
  222 + ***************************************************************************/
  223 +const GCTAOnOffObservation* GCTAOnOffObservations::at(const int& index) const
  224 +{
  225 + // Raise exception if index is out of range
  226 + if (index < 0 || index >= size()) {
  227 + throw GException::out_of_range(G_AT, index, 0, size()-1);
  228 + }
  229 +
  230 + // Return pointer
  231 + return m_obs[index];
  232 +}
  233 +
  234 +
  235 +/***********************************************************************//**
  236 + * @brief Set observation in container
  237 + *
  238 + * @param[in] index Observation index [0,...,size()-1].
  239 + * @param[in] obs Observation.
  240 + * @return Pointer to deep copy of observation.
  241 + *
  242 + * @exception GException::out_of_range
  243 + * Observation index is out of range.
  244 + * @exception GException::invalid_value
  245 + * Observation with same instrument and identifier already
  246 + * exists in container.
  247 + *
  248 + * Set a deep copy and observation @p obs at the specified @p index in the
  249 + * container.
  250 + ***************************************************************************/
  251 +GCTAOnOffObservation* GCTAOnOffObservations::set(const int& index, const GCTAOnOffObservation& obs)
  252 +{
  253 + // Compile option: raise exception if index is out of range
  254 + #if defined(G_RANGE_CHECK)
  255 + if (index < 0 || index >= size()) {
  256 + throw GException::out_of_range(G_SET, index, 0, size()-1);
  257 + }
  258 + #endif
  259 +
  260 + // Raise an exception if an observation with specified instrument and
  261 + // identifier already exists
  262 + int inx = get_index(obs.instrument(), obs.id());
  263 + if (inx != -1 && inx != index) {
  264 + std::string msg =
  265 + "Attempt to set \""+obs.instrument()+"\" observation with"
  266 + " identifier \""+obs.id()+"\" in observation container at"
  267 + " index "+gammalib::str(index)+", but an observation with the"
  268 + " same attributes exists already at index "+gammalib::str(inx)+
  269 + " in the container.\n"
  270 + "Every observation for a given instrument in the observation"
  271 + " container needs a unique identifier.";
  272 + throw GException::invalid_value(G_SET, msg);
  273 + }
  274 +
  275 + // Delete existing observation
  276 + if (m_obs[index] != NULL) delete m_obs[index];
  277 +
  278 + // Store pointer to a deep copy of the observation
  279 + m_obs[index] = obs.clone();
  280 +
  281 + // Return pointer
  282 + return m_obs[index];
  283 +}
  284 +
  285 +
  286 +/***********************************************************************//**
  287 + * @brief Append observation to container
  288 + *
  289 + * @param[in] obs Observation.
  290 + * @return Pointer to deep copy of observation.
  291 + *
  292 + * @exception GException::invalid_value
  293 + * Observation with same instrument and identifier already
  294 + * exists in container.
  295 + *
  296 + * Appends a deep copy of an observation to the container.
  297 + ***************************************************************************/
  298 +GCTAOnOffObservation* GCTAOnOffObservations::append(const GCTAOnOffObservation& obs)
  299 +{
  300 + // Raise an exception if an observation with specified instrument and
  301 + // identifier already exists
  302 + int inx = get_index(obs.instrument(), obs.id());
  303 + if (inx != -1) {
  304 + std::string msg =
  305 + "Attempt to append \""+obs.instrument()+"\" observation with"
  306 + " identifier \""+obs.id()+"\" to observation container, but an"
  307 + " observation with the same attributes exists already at"
  308 + " index "+gammalib::str(inx)+" in the container.\n"
  309 + "Every observation for a given instrument in the observation"
  310 + " container needs a unique identifier.";
  311 + throw GException::invalid_value(G_APPEND, msg);
  312 + }
  313 +
  314 + // Clone observation
  315 + GCTAOnOffObservation* ptr = obs.clone();
  316 +
  317 + // Append to list
  318 + m_obs.push_back(ptr);
  319 +
  320 + // Return pointer
  321 + return ptr;
  322 +}
  323 +
  324 +
  325 +/***********************************************************************//**
  326 + * @brief Insert observation into container
  327 + *
  328 + * @param[in] index Observation index [0,...,size()-1].
  329 + * @param[in] obs Observation.
  330 + * @return Pointer to deep copy of observation.
  331 + *
  332 + * @exception GException::out_of_range
  333 + * Observation index is out of range.
  334 + * @exception GException::invalid_value
  335 + * Observation with same instrument and identifier already
  336 + * exists in container.
  337 + *
  338 + * Inserts a deep copy of an observation into the container before the
  339 + * observation with the specified @p index.
  340 + ***************************************************************************/
  341 +GCTAOnOffObservation* GCTAOnOffObservations::insert(const int& index, const GCTAOnOffObservation& obs)
  342 +{
  343 + // Compile option: raise exception if index is out of range
  344 + #if defined(G_RANGE_CHECK)
  345 + if (isempty()) {
  346 + if (index > 0) {
  347 + throw GException::out_of_range(G_INSERT, index, 0, size()-1);
  348 + }
  349 + }
  350 + else {
  351 + if (index < 0 || index >= size()) {
  352 + throw GException::out_of_range(G_INSERT, index, 0, size()-1);
  353 + }
  354 + }
  355 + #endif
  356 +
  357 + // Raise an exception if an observation with specified instrument and
  358 + // identifier already exists
  359 + int inx = get_index(obs.instrument(), obs.id());
  360 + if (inx != -1 && inx != index) {
  361 + std::string msg =
  362 + "Attempt to insert \""+obs.instrument()+"\" observation with"
  363 + " identifier \""+obs.id()+"\" in observation container before"
  364 + " index "+gammalib::str(index)+", but an observation with the"
  365 + " same attributes exists already at index "+gammalib::str(inx)+
  366 + " in the container.\n"
  367 + "Every observation for a given instrument in the observation"
  368 + " container needs a unique identifier.";
  369 + throw GException::invalid_value(G_INSERT, msg);
  370 + }
  371 +
  372 + // Clone observation
  373 + GCTAOnOffObservation* ptr = obs.clone();
  374 +
  375 + // Clone observation and insert into list
  376 + m_obs.insert(m_obs.begin()+index, ptr);
  377 +
  378 + // Return pointer
  379 + return ptr;
  380 +}
  381 +
  382 +
  383 +/***********************************************************************//**
  384 + * @brief Remove observation from container
  385 + *
  386 + * @param[in] index Observation index [0,...,size()-1].
  387 + *
  388 + * @exception GException::out_of_range
  389 + * Observation index is out of range.
  390 + *
  391 + * Removes observation of specified @p index from the container.
  392 + ***************************************************************************/
  393 +void GCTAOnOffObservations::remove(const int& index)
  394 +{
  395 + // Compile option: If index is outside boundary then raise exception
  396 + #if defined(G_RANGE_CHECK)
  397 + if (index < 0 || index >= size()) {
  398 + throw GException::out_of_range(G_REMOVE, index, 0, size()-1);
  399 + }
  400 + #endif
  401 +
  402 + // Erase observation component from container
  403 + m_obs.erase(m_obs.begin() + index);
  404 +
  405 + // Return
  406 + return;
  407 +}
  408 +
  409 +
  410 +/***********************************************************************//**
  411 + * @brief Append observations from observation container
  412 + *
  413 + * @param[in] obs Observations.
  414 + *
  415 + * @exception GException::invalid_value
  416 + * Observation with same instrument and identifier already
  417 + * exists in container.
  418 + *
  419 + * Appends deep copies of observations to the container.
  420 + ***************************************************************************/
  421 +void GCTAOnOffObservations::extend(const GCTAOnOffObservations& obs)
  422 +{
  423 + // Do nothing if observation container is empty
  424 + if (!obs.isempty()) {
  425 +
  426 + // Get size. Note that we extract the size first to avoid an
  427 + // endless loop that arises when a container is appended to
  428 + // itself.
  429 + int num = obs.size();
  430 +
  431 + // Reserve enough space
  432 + reserve(size() + num);
  433 +
  434 + // Loop over all observations, clone them and append them to the
  435 + // list
  436 + for (int i = 0; i < num; ++i) {
  437 +
  438 + // Raise an exception if an observation with specified
  439 + // instrument and identifier already exists
  440 + int inx = get_index(obs[i]->instrument(), obs[i]->id());
  441 + if (inx != -1) {
  442 + std::string msg =
  443 + "Attempt to append \""+obs[i]->instrument()+"\""
  444 + " observation with identifier \""+obs[i]->id()+"\""
  445 + " to observation container, but an observation with the"
  446 + " same attributes exists already at index "+
  447 + gammalib::str(inx)+" in the container.\n"
  448 + "Every observation for a given instrument in the"
  449 + " observation container needs a unique identifier.";
  450 + throw GException::invalid_value(G_EXTEND, msg);
  451 + }
  452 +
  453 + // Append observation to container
  454 + m_obs.push_back(obs[i]->clone());
  455 + }
  456 +
  457 + } // endif: observation container was not empty
  458 +
  459 + // Return
  460 + return;
  461 +}
  462 +
  463 +
  464 +/***********************************************************************//**
  465 + * @brief Signals if observation exists
  466 + *
  467 + * @param[in] instrument Instrument.
  468 + * @param[in] id Observation identifier.
  469 + * @return True if observation with specified @p instrument and identifier
  470 + * @p id exists.
  471 + *
  472 + * Searches all observations for a match with the specified @p instrument
  473 + * and @p identifier. If the specified attributes have been found, true is
  474 + * returned.
  475 + ***************************************************************************/
  476 +bool GCTAOnOffObservations::contains(const std::string& instrument,
  477 + const std::string& id) const
  478 +{
  479 + // Get observation index
  480 + int index = get_index(instrument, id);
  481 +
  482 + // Return
  483 + return (index != -1);
  484 +}
  485 +
  486 +
  487 +/***********************************************************************//**
  488 + * @brief Load observations from XML file
  489 + *
  490 + * @param[in] filename XML filename.
  491 + *
  492 + * Loads observation from a XML file into the container. Please refer to the
  493 + * read() method for more information about the structure of the XML file.
  494 + ***************************************************************************/
  495 +void GCTAOnOffObservations::load(const std::string& filename)
  496 +{
  497 + // Clear any existing observations
  498 + clear();
  499 +
  500 + // Load XML document
  501 + GXml xml(filename);
  502 +
  503 + // Read observations from XML document
  504 + read(xml);
  505 +
  506 + // Return
  507 + return;
  508 +}
  509 +
  510 +
  511 +/***********************************************************************//**
  512 + * @brief Save observations into XML file
  513 + *
  514 + * @param[in] filename XML filename.
  515 + *
  516 + * Saves observations into a XML file. Please refer to the read() method for
  517 + * more information about the structure of the XML file.
  518 + ***************************************************************************/
  519 +void GCTAOnOffObservations::save(const std::string& filename) const
  520 +{
  521 + // Declare empty XML document
  522 + GXml xml;
  523 +
  524 + // Write observations into XML file
  525 + write(xml);
  526 +
  527 + // Save XML document
  528 + xml.save(filename);
  529 +
  530 + // Return
  531 + return;
  532 +}
  533 +
  534 +
  535 +/***********************************************************************//**
  536 + * @brief Read observations from XML document
  537 + *
  538 + * @param[in] xml XML document.
  539 + *
  540 + * @exception GException::invalid_instrument
  541 + * Invalid instrument encountered in XML file.
  542 + *
  543 + * Reads observations from the first observation list that is found in the
  544 + * XML document. The decoding of the instrument specific observation
  545 + * definition is done within the observation's GCTAOnOffObservation::read() method.
  546 + * The following file structure is expected:
  547 + *
  548 + * <observation_list title="observation library">
  549 + * <observation name="..." id="..." instrument="...">
  550 + * ...
  551 + * </observation>
  552 + * <observation name="..." id="..." instrument="...">
  553 + * ...
  554 + * </observation>
  555 + * ...
  556 + * </observation_list>
  557 + *
  558 + * The @p name and @p id attributes allow for a unique identification of an
  559 + * observation within the observation container. The @p instrument
  560 + * attributes specifies the instrument to which the observation applies.
  561 + * At the moment there is only CTA but more instruments may be implemented
  562 + * later and require the need for a registry (see GObservations.cpp).
  563 + *
  564 + * The structure within the @p observation tag is defined by the instrument
  565 + * specific GCTAOnOffObservation class.
  566 + *
  567 + ***************************************************************************/
  568 +void GCTAOnOffObservations::read(const GXml& xml)
  569 +{
  570 + // Get pointer on observation library
  571 + const GXmlElement* lib = xml.element("observation_list", 0);
  572 +
  573 + // Loop over all observations
  574 + int n = lib->elements("observation");
  575 + for (int i = 0; i < n; ++i) {
  576 +
  577 + // Get pointer on observation
  578 + const GXmlElement* obs = lib->element("observation", i);
  579 +
  580 + // Get attributes
  581 + std::string name = obs->attribute("name");
  582 + std::string id = obs->attribute("id");
  583 + std::string instrument = obs->attribute("instrument");
  584 +
  585 + // Allocate observation
  586 + // (only CTA at the moment, implement registry if more)
  587 + GCTAOnOffObservation* ptr = new GCTAOnOffObservation;
  588 +
  589 + // If observation is valid then read its definition from XML file
  590 + if (ptr != NULL) {
  591 +
  592 + // Read definition
  593 + ptr->read(*obs);
  594 +
  595 + // Set attributes
  596 + ptr->name(name);
  597 + ptr->id(id);
  598 +
  599 + } // endif: observation was valid
  600 +
  601 + // ... otherwise throw an exception
  602 + else {
  603 + std::string msg =
  604 + "Problem allocating GCTAOnOffObservation object.";
  605 + throw GException::invalid_instrument(G_READ, msg);
  606 + }
  607 +
  608 + // Append observation to container
  609 + append(*ptr);
  610 +
  611 + // Free observation (the append method made a deep copy)
  612 + delete ptr;
  613 +
  614 + } // endfor: looped over all observations
  615 +
  616 + // Return
  617 + return;
  618 +}
  619 +
  620 +
  621 +/***********************************************************************//**
  622 + * @brief Write observations into XML document
  623 + *
  624 + * @param[in] xml XML document.
  625 + *
  626 + * Write observations into the first observation library that is found in the
  627 + * XML document. In case that no observation library exists, one is added to
  628 + * the document. Please refer to the read() method for more information
  629 + * about the structure of the XML document.
  630 + ***************************************************************************/
  631 +void GCTAOnOffObservations::write(GXml& xml) const
  632 +{
  633 + // If there is no observation library then append one
  634 + if (xml.elements("observation_list") == 0) {
  635 + xml.append(GXmlElement("observation_list title=\"observation list\""));
  636 + }
  637 +
  638 + // Get pointer on observation library
  639 + GXmlElement* lib = xml.element("observation_list", 0);
  640 +
  641 + // Write all observations into library
  642 + for (int i = 0; i < size(); ++i) {
  643 +
  644 + // Initialise pointer on observation
  645 + GXmlElement* obs = NULL;
  646 +
  647 + // Search corresponding observation
  648 + int n = xml.elements("observation");
  649 + for (int k = 0; k < n; ++k) {
  650 + GXmlElement* element = xml.element("observation", k);
  651 + if (element->attribute("name") == m_obs[i]->name() &&
  652 + element->attribute("id") == m_obs[i]->id() &&
  653 + element->attribute("instrument") == m_obs[i]->instrument()) {
  654 + obs = element;
  655 + break;
  656 + }
  657 + }
  658 +
  659 + // If no observation with corresponding name, ID and instrument was
  660 + // found then append one now
  661 + if (obs == NULL) {
  662 + obs = lib->append("observation");
  663 + obs->attribute("name", m_obs[i]->name());
  664 + obs->attribute("id", m_obs[i]->id());
  665 + obs->attribute("instrument", m_obs[i]->instrument());
  666 + }
  667 +
  668 + // Write now observation
  669 + m_obs[i]->write(*obs);
  670 +
  671 + } // endfor: looped over all observaitons
  672 +
  673 + // Return
  674 + return;
  675 +}
  676 +
  677 +
  678 +/***********************************************************************//**
  679 + * @brief Load models from XML file
  680 + *
  681 + * @param[in] filename XML filename.
  682 + *
  683 + * Loads the models from a XML file. Please refer to the GModels::read()
  684 + * method for more information about the expected structure of the XML
  685 + * file.
  686 + ***************************************************************************/
  687 +void GCTAOnOffObservations::models(const std::string& filename)
  688 +{
  689 + // Load models
  690 + m_models.load(filename);
  691 +
  692 + // Return
  693 + return;
  694 +}
  695 +
  696 +
  697 +/***********************************************************************//**
  698 + * @brief Print observation list information
  699 + *
  700 + * @param[in] chatter Chattiness (defaults to NORMAL).
  701 + * @return String containing observation list information
  702 + ***************************************************************************/
  703 +std::string GCTAOnOffObservations::print(const GChatter& chatter) const
  704 +{
  705 + // Initialise result string
  706 + std::string result;
  707 +
  708 + // Continue only if chatter is not silent
  709 + if (chatter != SILENT) {
  710 +
  711 + // Append header
  712 + result.append("=== GCTAOnOffObservations ===");
  713 +
  714 + // Append information
  715 + result.append("\n"+gammalib::parformat("Number of observations"));
  716 + result.append(gammalib::str(size()));
  717 +
  718 + // NORMAL: Append observations
  719 + if (chatter >= NORMAL) {
  720 + for (int i = 0; i < size(); ++i) {
  721 + result.append("\n");
  722 + result.append((*this)[i]->print());
  723 + }
  724 + } // endif: chatter was normal
  725 +
  726 + // EXPLICIT: Append models
  727 + if (chatter >= EXPLICIT) {
  728 + result.append("\n"+m_models.print());
  729 + } // endif: chatter was explicit
  730 +
  731 + } // endif: chatter was not silent
  732 +
  733 + // Return result
  734 + return result;
  735 +}
  736 +
  737 +
  738 +/*==========================================================================
  739 + = =
  740 + = Private methods =
  741 + = =
  742 + ==========================================================================*/
  743 +
  744 +/***********************************************************************//**
  745 + * @brief Initialise class members
  746 + ***************************************************************************/
  747 +void GCTAOnOffObservations::init_members(void)
  748 +{
  749 + // Initialise members
  750 + m_obs.clear();
  751 + m_models.clear();
  752 +
  753 + // Return
  754 + return;
  755 +}
  756 +
  757 +
  758 +/***********************************************************************//**
  759 + * @brief Copy class members
  760 + *
  761 + * @param[in] obs Observation container.
  762 + *
  763 + * Copies all class members. Deep copies are created for the observations,
  764 + * which is what is needed to have a real physical copy of the members.
  765 + ***************************************************************************/
  766 +void GCTAOnOffObservations::copy_members(const GCTAOnOffObservations& obs)
  767 +{
  768 + // Copy attributes
  769 + m_models = obs.m_models;
  770 +
  771 + // Copy observations
  772 + m_obs.clear();
  773 + for (int i = 0; i < obs.m_obs.size(); ++i) {
  774 + m_obs.push_back((obs.m_obs[i]->clone()));
  775 + }
  776 +
  777 + // Return
  778 + return;
  779 +}
  780 +
  781 +
  782 +/***********************************************************************//**
  783 + * @brief Delete class members
  784 + *
  785 + * Deallocates all observations. Since container classes that hold pointers
  786 + * need to handle the proper deallocation of memory, we loop here over all
  787 + * pointers and make sure that we deallocate all observations in the
  788 + * container.
  789 + ***************************************************************************/
  790 +void GCTAOnOffObservations::free_members(void)
  791 +{
  792 + // Free observations
  793 + for (int i = 0; i < m_obs.size(); ++i) {
  794 + delete m_obs[i];
  795 + m_obs[i] = NULL;
  796 + }
  797 +
  798 + // Return
  799 + return;
  800 +}
  801 +
  802 +
  803 +/***********************************************************************//**
  804 + * @brief Return observation index by instrument and identifier
  805 + *
  806 + * @param[in] instrument Instrument.
  807 + * @param[in] id Observation identifier.
  808 + * @return Observation index (-1 of instrument and observation identifier
  809 + * has not been found)
  810 + *
  811 + * Returns observation index based on the specified @p instrument and
  812 + * observation identifier @p id. If no observation with the specified
  813 + * attributes has been found, the method returns -1.
  814 + ***************************************************************************/
  815 +int GCTAOnOffObservations::get_index(const std::string& instrument,
  816 + const std::string& id) const
  817 +{
  818 + // Initialise index
  819 + int index = -1;
  820 +
  821 + // Search observation with specified instrument and id
  822 + for (int i = 0; i < size(); ++i) {
  823 + if ((m_obs[i]->instrument() == instrument) &&
  824 + (m_obs[i]->id() == id)) {
  825 + index = i;
  826 + break;
  827 + }
  828 + }
  829 +
  830 + // Return index
  831 + return index;
  832 +}
... ...