Simon Bonnefoy / gammalib

Commit 513f99ba18ae7fd4685bfaa374da3a1181ff816f

Authored by Jürgen Knödlseder
1 parent 642d66bb

Implement skeleton for general CTA background model (#2694)

Inline Side-by-side
Showing 16 changed files with 2373 additions and 43 deletions
inst/cta/Makefile.am
  View file @513f99b
... ... @@ -74,6 +74,9 @@ sources = src/GCTASupport.cpp \
74 74 src/GCTAInstDir.cpp \
75 75 src/GCTARoi.cpp \
76 76 src/GCTAPointing.cpp \
  77 + src/GCTAModelBackground.cpp \
  78 + src/GCTAModelSpatial.cpp \
  79 + src/GCTAModelSpatialRegistry.cpp \
77 80 src/GCTAModelCubeBackground.cpp \
78 81 src/GCTAModelIrfBackground.cpp \
79 82 src/GCTAModelAeffBackground.cpp \
... ... @@ -124,6 +127,9 @@ pkginclude_HEADERS = include/GCTATypemaps.hpp \
124 127 include/GCTACubeSource.hpp \
125 128 include/GCTACubeSourcePoint.hpp \
126 129 include/GCTACubeSourceDiffuse.hpp \
  130 + include/GCTAModelBackground.hpp \
  131 + include/GCTAModelSpatial.hpp \
  132 + include/GCTAModelSpatialRegistry.hpp \
127 133 include/GCTAModelCubeBackground.hpp \
128 134 include/GCTAModelIrfBackground.hpp \
129 135 include/GCTAModelAeffBackground.hpp \
... ...
inst/cta/include/GCTALib.hpp
  View file @513f99b
1 1 /***************************************************************************
2 2 * GCTALib.hpp - CTA Support Header files *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -72,6 +72,9 @@
72 72 #include "GCTACubeSource.hpp"
73 73 #include "GCTACubeSourcePoint.hpp"
74 74 #include "GCTACubeSourceDiffuse.hpp"
  75 +#include "GCTAModelBackground.hpp"
  76 +#include "GCTAModelSpatial.hpp"
  77 +#include "GCTAModelSpatialRegistry.hpp"
75 78 #include "GCTAModelCubeBackground.hpp"
76 79 #include "GCTAModelIrfBackground.hpp"
77 80 #include "GCTAModelAeffBackground.hpp"
... ...
inst/cta/include/GCTAModelBackground.hpp 0 → 100644
  View file @513f99b
  1 +/***************************************************************************
  2 + * GCTAModelBackground.hpp - Background model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAModelBackground.hpp
  23 + * @brief Background model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCTAMODELBACKGROUND_HPP
  28 +#define GCTAMODELBACKGROUND_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include <cmath>
  33 +#include "GModelData.hpp"
  34 +#include "GModelSpectral.hpp"
  35 +#include "GModelTemporal.hpp"
  36 +#include "GFunction.hpp"
  37 +#include "GCTAModelSpatial.hpp"
  38 +#include "GCTAEventList.hpp"
  39 +
  40 +/* __ Forward declarations _______________________________________________ */
  41 +class GEvent;
  42 +class GObservation;
  43 +class GModelPar;
  44 +class GXmlElement;
  45 +
  46 +
  47 +/***********************************************************************//**
  48 + * @class GCTAModelBackground
  49 + *
  50 + * @brief Background model class
  51 + *
  52 + * This class implements a background model for CTA.
  53 + ***************************************************************************/
  54 +class GCTAModelBackground : public GModelData {
  55 +
  56 +public:
  57 + // Constructors and destructors
  58 + GCTAModelBackground(void);
  59 + explicit GCTAModelBackground(const GXmlElement& xml);
  60 + GCTAModelBackground(const GCTAModelSpatial& spatial,
  61 + const GModelSpectral& spectral);
  62 + GCTAModelBackground(const GCTAModelSpatial& spatial,
  63 + const GModelSpectral& spectral,
  64 + const GModelTemporal& temporal);
  65 + GCTAModelBackground(const GCTAModelBackground& model);
  66 + virtual ~GCTAModelBackground(void);
  67 +
  68 + // Operators
  69 + virtual GCTAModelBackground& operator=(const GCTAModelBackground& model);
  70 +
  71 + // Implemented pure virtual methods
  72 + virtual void clear(void);
  73 + virtual GCTAModelBackground* clone(void) const;
  74 + virtual std::string classname(void) const;
  75 + virtual std::string type(void) const;
  76 + virtual bool is_constant(void) const;
  77 + virtual double eval(const GEvent& event,
  78 + const GObservation& obs,
  79 + const bool& gradients = false) const;
  80 + virtual double npred(const GEnergy& obsEng,
  81 + const GTime& obsTime,
  82 + const GObservation& obs) const;
  83 + virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
  84 + virtual void read(const GXmlElement& xml);
  85 + virtual void write(GXmlElement& xml) const;
  86 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  87 +
  88 + // Other methods
  89 + GCTAModelSpatial* spatial(void) const;
  90 + GModelSpectral* spectral(void) const;
  91 + GModelTemporal* temporal(void) const;
  92 + void spatial(const GCTAModelSpatial* spatial);
  93 + void spectral(const GModelSpectral* spectral);
  94 + void temporal(const GModelTemporal* temporal);
  95 +
  96 +protected:
  97 + // Protected methods
  98 + void init_members(void);
  99 + void copy_members(const GCTAModelBackground& model);
  100 + void free_members(void);
  101 + void set_pointers(void);
  102 + GCTAModelSpatial* xml_spatial(const GXmlElement& spatial) const;
  103 + GModelSpectral* xml_spectral(const GXmlElement& spectral) const;
  104 + GModelTemporal* xml_temporal(const GXmlElement& temporal) const;
  105 + bool valid_model(void) const;
  106 +
  107 + // ROI integration kernel
  108 + class roi_kern : public GFunction {
  109 + public:
  110 + roi_kern(const GCTAModelSpatial* parent, const double& roi, const double& dist) :
  111 + m_parent(parent),
  112 + m_roi(roi),
  113 + m_cosroi(std::cos(roi)),
  114 + m_dist(dist),
  115 + m_cosdist(std::cos(dist)),
  116 + m_sindist(std::sin(dist)) { }
  117 + double eval(const double& r);
  118 + protected:
  119 + const GCTAModelSpatial* m_parent; //!< Pointer to radial model
  120 + double m_roi; //!< ROI radius in radians
  121 + double m_cosroi; //!< Cosine of ROI radius
  122 + double m_dist; //!< Distance between pointing and ROI centre in radians
  123 + double m_cosdist; //!< Cosine of distance
  124 + double m_sindist; //!< Sinus of distance
  125 + };
  126 +
  127 + // Proteced data members
  128 + GCTAModelSpatial* m_spatial; //!< Spatial model
  129 + GModelSpectral* m_spectral; //!< Spectral model
  130 + GModelTemporal* m_temporal; //!< Temporal model
  131 +};
  132 +
  133 +
  134 +/***********************************************************************//**
  135 + * @brief Return class name
  136 + *
  137 + * @return String containing the class name ("GCTAModelBackground").
  138 + ***************************************************************************/
  139 +inline
  140 +std::string GCTAModelBackground::classname(void) const
  141 +{
  142 + return ("GCTAModelBackground");
  143 +}
  144 +
  145 +
  146 +/***********************************************************************//**
  147 + * @brief Return model type
  148 + *
  149 + * @return Model type.
  150 + *
  151 + * Returns "CTABackground" as the type of the model.
  152 + ***************************************************************************/
  153 +inline
  154 +std::string GCTAModelBackground::type(void) const
  155 +{
  156 + return ("CTABackground");
  157 +}
  158 +
  159 +
  160 +/***********************************************************************//**
  161 + * @brief Signals if model is temporally constant
  162 + *
  163 + * @return True if model is temporally constant, false otherwise.
  164 + *
  165 + * Signals if the model is temporally constant. A temporally constant model
  166 + * is a model that has a temporal component of type "Constant".
  167 + ***************************************************************************/
  168 +inline
  169 +bool GCTAModelBackground::is_constant(void) const
  170 +{
  171 + return (m_temporal != NULL && m_temporal->type() == "Constant");
  172 +}
  173 +
  174 +
  175 +/***********************************************************************//**
  176 + * @brief Return spatial model component
  177 + *
  178 + * @return Pointer to spatial model component.
  179 + *
  180 + * Returns a pointer to the spatial model component of the model. The pointer
  181 + * is of type GCTAModelSpatial. Note that a NULL pointer may be returned if
  182 + * the model has no spatial model component.
  183 + ***************************************************************************/
  184 +inline
  185 +GCTAModelSpatial* GCTAModelBackground::spatial(void) const
  186 +{
  187 + return (m_spatial);
  188 +}
  189 +
  190 +
  191 +/***********************************************************************//**
  192 + * @brief Return spectral model component
  193 + *
  194 + * @return Pointer to spectral model component.
  195 + *
  196 + * Returns a pointer to the spectral model component of the model. The
  197 + * pointer is of type GModelSpectral. Note that a NULL pointer may be
  198 + * returned if the model has no spectral model component.
  199 + ***************************************************************************/
  200 +inline
  201 +GModelSpectral* GCTAModelBackground::spectral(void) const
  202 +{
  203 + return (m_spectral);
  204 +}
  205 +
  206 +
  207 +/***********************************************************************//**
  208 + * @brief Return temporal model component
  209 + *
  210 + * @return Pointer to temporal model component.
  211 + *
  212 + * Returns a pointer to the temporal model component of the model. The
  213 + * pointer is of type GModelTemporal. Note that a NULL pointer may be
  214 + * returned if the model has no temporal model component.
  215 + ***************************************************************************/
  216 +inline
  217 +GModelTemporal* GCTAModelBackground::temporal(void) const
  218 +{
  219 + return (m_temporal);
  220 +}
  221 +
  222 +#endif /* GCTAMODELBACKGROUND_HPP */
... ...
inst/cta/include/GCTAModelRadial.hpp
  View file @513f99b
... ... @@ -29,11 +29,13 @@
29 29  
30 30 /* __ Includes ___________________________________________________________ */
31 31 #include <string>
32   -#include "GBase.hpp"
  32 +#include "GCTAModelSpatial.hpp"
33 33 #include "GModelPar.hpp"
34   -#include "GXmlElement.hpp"
35   -#include "GRan.hpp"
36   -#include "GCTAInstDir.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GRan;
  37 +class GXmlElement;
  38 +class GCTAInstDir;
37 39  
38 40  
39 41 /***********************************************************************//**
... ... @@ -44,7 +46,7 @@
44 46 * This class implements the radial component of the CTA radial acceptance
45 47 * model.
46 48 ***************************************************************************/
47   -class GCTAModelRadial : public GBase {
  49 +class GCTAModelRadial : public GCTAModelSpatial {
48 50  
49 51 public:
50 52 // Constructors and destructors
... ... @@ -59,19 +61,30 @@ public:
59 61 virtual GModelPar& operator[](const std::string& name);
60 62 virtual const GModelPar& operator[](const std::string& name) const;
61 63  
  64 + // Implemented virtual methods
  65 + virtual double eval(const GCTAInstDir& dir,
  66 + const GEnergy& energy,
  67 + const GTime& time,
  68 + const bool& gradients = false) const;
  69 + virtual GCTAInstDir mc(const GEnergy& energy,
  70 + const GTime& time,
  71 + GRan& ran) const;
  72 +
62 73 // Pure virtual methods
63 74 virtual void clear(void) = 0;
64 75 virtual GCTAModelRadial* clone(void) const = 0;
65 76 virtual std::string classname(void) const = 0;
66 77 virtual std::string type(void) const = 0;
67   - virtual double eval(const double& offset,
68   - const bool& gradients = false) const = 0;
69   - virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0;
70 78 virtual double omega(void) const = 0;
71 79 virtual void read(const GXmlElement& xml) = 0;
72 80 virtual void write(GXmlElement& xml) const = 0;
73 81 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
74 82  
  83 + // Derived class pure virtual methods
  84 + virtual double eval(const double& offset,
  85 + const bool& gradients = false) const = 0;
  86 + virtual GCTAInstDir mc(const GCTAInstDir& dir, GRan& ran) const = 0;
  87 +
75 88 // Methods
76 89 int size(void) const;
77 90  
... ...
inst/cta/include/GCTAModelSpatial.hpp 0 → 100644
  View file @513f99b
  1 +/***************************************************************************
  2 + * GCTAModelSpatial.hpp - Spatial model abstract base class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAModelSpatial.hpp
  23 + * @brief Abstract spatial model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCTAMODELSPATIAL_HPP
  28 +#define GCTAMODELSPATIAL_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GBase.hpp"
  33 +#include "GModelPar.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GRan;
  37 +class GTime;
  38 +class GXmlElement;
  39 +class GCTAInstDir;
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class GCTAModelSpatial
  44 + *
  45 + * @brief Abstract spatial model class
  46 + *
  47 + * This class implements the spatial component of the CTA background model.
  48 + ***************************************************************************/
  49 +class GCTAModelSpatial : public GBase {
  50 +
  51 +public:
  52 + // Constructors and destructors
  53 + GCTAModelSpatial(void);
  54 + GCTAModelSpatial(const GCTAModelSpatial& model);
  55 + virtual ~GCTAModelSpatial(void);
  56 +
  57 + // Operators
  58 + virtual GCTAModelSpatial& operator=(const GCTAModelSpatial& model);
  59 + virtual GModelPar& operator[](const int& index);
  60 + virtual const GModelPar& operator[](const int& index) const;
  61 + virtual GModelPar& operator[](const std::string& name);
  62 + virtual const GModelPar& operator[](const std::string& name) const;
  63 +
  64 + // Pure virtual methods
  65 + virtual void clear(void) = 0;
  66 + virtual GCTAModelSpatial* clone(void) const = 0;
  67 + virtual std::string classname(void) const = 0;
  68 + virtual std::string type(void) const = 0;
  69 + virtual double eval(const GCTAInstDir& dir,
  70 + const GEnergy& energy,
  71 + const GTime& time,
  72 + const bool& gradients = false) const = 0;
  73 + virtual GCTAInstDir mc(const GEnergy& energy,
  74 + const GTime& time,
  75 + GRan& ran) const = 0;
  76 + virtual double omega(void) const = 0;
  77 + virtual void read(const GXmlElement& xml) = 0;
  78 + virtual void write(GXmlElement& xml) const = 0;
  79 + virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
  80 +
  81 + // Methods
  82 + int size(void) const;
  83 +
  84 +protected:
  85 + // Protected methods
  86 + void init_members(void);
  87 + void copy_members(const GCTAModelSpatial& model);
  88 + void free_members(void);
  89 +
  90 + // Proteced members
  91 + std::vector<GModelPar*> m_pars; //!< Parameter pointers
  92 +};
  93 +
  94 +
  95 +/***********************************************************************//**
  96 + * @brief Return number of model parameters
  97 + *
  98 + * @return Number of model parameters.
  99 + ***************************************************************************/
  100 +inline
  101 +int GCTAModelSpatial::size(void) const
  102 +{
  103 + return ((int)m_pars.size());
  104 +}
  105 +
  106 +#endif /* GCTAMODELSPATIAL_HPP */
... ...
inst/cta/include/GCTAModelSpatialRegistry.hpp 0 → 100644
  View file @513f99b
  1 +/***************************************************************************
  2 + * GCTAModelSpatialRegistry.hpp - Spatial model registry class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAModelSpatialRegistry.hpp
  23 + * @brief Spatial model registry class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCTAMODELSPATIALREGISTRY_HPP
  28 +#define GCTAMODELSPATIALREGISTRY_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GRegistry.hpp"
  33 +
  34 +/* __ Forward declarations _______________________________________________ */
  35 +class GXmlElement;
  36 +class GCTAModelSpatial;
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class GCTAModelSpatialRegistry
  41 + *
  42 + * @brief Interface definition for the spatial model registry class
  43 + *
  44 + * The registry class allows the registration of spatial models that are not
  45 + * necessarily compiled into the GammaLib. It uses the static members
  46 + * m_number and m_models which are allocated globally to keep track of
  47 + * spatial models that are available throughout all linked libraries. To
  48 + * register a spatial model it is sufficient to add
  49 + *
  50 + * const GCTAModelSpatialXXX g_spatial_XXX_seed;
  51 + * const GCTAModelSpatialRegistry g_spatial_XXX_registry(&g_spatial_XXX_seed);
  52 + *
  53 + * at the top of the .cpp file of the spatial model. Here, XXX is a unique
  54 + * name that describes the model.
  55 + ***************************************************************************/
  56 +class GCTAModelSpatialRegistry : public GRegistry {
  57 +
  58 +public:
  59 + // Constructors and destructors
  60 + GCTAModelSpatialRegistry(void);
  61 + GCTAModelSpatialRegistry(const GCTAModelSpatial* model);
  62 + GCTAModelSpatialRegistry(const GCTAModelSpatialRegistry& registry);
  63 + virtual ~GCTAModelSpatialRegistry(void);
  64 +
  65 + // Operators
  66 + GCTAModelSpatialRegistry& operator=(const GCTAModelSpatialRegistry& registry);
  67 +
  68 + // Methods
  69 + std::string classname(void) const;
  70 + int size(void) const;
  71 + GCTAModelSpatial* alloc(const GXmlElement& xml) const;
  72 + std::string name(const int& index) const;
  73 + std::string print(const GChatter& chatter = NORMAL) const;
  74 +
  75 +protected:
  76 + // Protected methods
  77 + void init_members(void);
  78 + void copy_members(const GCTAModelSpatialRegistry& registry);
  79 + void free_members(void);
  80 +
  81 +private:
  82 + // Private members (the private members have been implement as static
  83 + // methods to avoid the static initialization order fiasco of static
  84 + // members; using static methods we follow the "construct on first use
  85 + // idiom")
  86 + // Number of models in registry
  87 + static int& number() {
  88 + static int m_number = 0;
  89 + return m_number;
  90 + };
  91 + // Pointer to seed models
  92 + static GRegistryPointer<const GCTAModelSpatial*>& models() {
  93 + static GRegistryPointer<const GCTAModelSpatial*> m_models;
  94 + return m_models;
  95 + };
  96 +};
  97 +
  98 +
  99 +/***********************************************************************//**
  100 + * @brief Return class name
  101 + *
  102 + * @return String containing the class name ("GCTAModelSpatialRegistry").
  103 + ***************************************************************************/
  104 +inline
  105 +std::string GCTAModelSpatialRegistry::classname(void) const
  106 +{
  107 + return ("GCTAModelSpatialRegistry");
  108 +}
  109 +
  110 +
  111 +/***********************************************************************//**
  112 + * @brief Return number of registered models
  113 + *
  114 + * @return Number of registered models.
  115 + *
  116 + * Returns the number of registered model.
  117 + ***************************************************************************/
  118 +inline
  119 +int GCTAModelSpatialRegistry::size(void) const
  120 +{
  121 + return number();
  122 +}
  123 +
  124 +#endif /* GCTAMODELSPATIALREGISTRY_HPP */
... ...
inst/cta/src/GCTAModelBackground.cpp 0 → 100644
  View file @513f99b
  1 +/***************************************************************************
  2 + * GCTAModelBackground.cpp - Background model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAModelBackground.cpp
  23 + * @brief Background model class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GTools.hpp"
  33 +#include "GMath.hpp"
  34 +#include "GModelRegistry.hpp"
  35 +#include "GModelSpectralRegistry.hpp"
  36 +#include "GModelTemporalRegistry.hpp"
  37 +#include "GModelTemporalConst.hpp"
  38 +#include "GIntegral.hpp"
  39 +#include "GCTAModelSpatialRegistry.hpp"
  40 +#include "GCTAModelBackground.hpp"
  41 +#include "GCTAObservation.hpp"
  42 +#include "GCTAPointing.hpp"
  43 +#include "GCTAInstDir.hpp"
  44 +#include "GCTARoi.hpp"
  45 +#include "GCTAException.hpp"
  46 +#include "GCTASupport.hpp"
  47 +
  48 +/* __ Constants __________________________________________________________ */
  49 +
  50 +/* __ Globals ____________________________________________________________ */
  51 +const GCTAModelBackground g_cta_background_seed;
  52 +const GModelRegistry g_cta_background_registry(&g_cta_background_seed);
  53 +
  54 +/* __ Method name definitions ____________________________________________ */
  55 +#define G_EVAL "GCTAModelBackground::eval(GEvent&, GObservation&, bool&)"
  56 +#define G_NPRED "GCTAModelBackground::npred(GEnergy&, GTime&, GObservation&)"
  57 +#define G_MC "GCTAModelBackground::mc(GObservation&, GRan&)"
  58 +#define G_XML_SPATIAL "GCTAModelBackground::xml_spatial(GXmlElement&)"
  59 +#define G_XML_SPECTRAL "GCTAModelBackground::xml_spectral(GXmlElement&)"
  60 +#define G_XML_TEMPORAL "GCTAModelBackground::xml_temporal(GXmlElement&)"
  61 +
  62 +/* __ Macros _____________________________________________________________ */
  63 +
  64 +/* __ Coding definitions _________________________________________________ */
  65 +
  66 +/* __ Debug definitions __________________________________________________ */
  67 +//#define G_DUMP_MC //!< Dump MC information
  68 +
  69 +
  70 +/*==========================================================================
  71 + = =
  72 + = Constructors/destructors =
  73 + = =
  74 + ==========================================================================*/
  75 +
  76 +/***********************************************************************//**
  77 + * @brief Void constructor
  78 + *
  79 + * Constructs an empty background model.
  80 + ***************************************************************************/
  81 +GCTAModelBackground::GCTAModelBackground(void) : GModelData()
  82 +{
  83 + // Initialise members
  84 + init_members();
  85 +
  86 + // Return
  87 + return;
  88 +}
  89 +
  90 +
  91 +/***********************************************************************//**
  92 + * @brief XML constructor
  93 + *
  94 + * @param[in] xml XML element.
  95 + *
  96 + * Constructs a background model from the information that is found in a XML
  97 + * element. Please refer to the read() method to learn about the expected
  98 + * structure of the XML element.
  99 + ***************************************************************************/
  100 +GCTAModelBackground::GCTAModelBackground(const GXmlElement& xml) :
  101 + GModelData(xml)
  102 +{
  103 + // Initialise members
  104 + init_members();
  105 +
  106 + // Read XML
  107 + read(xml);
  108 +
  109 + // Set parameter pointers
  110 + set_pointers();
  111 +
  112 + // Return
  113 + return;
  114 +}
  115 +
  116 +
  117 +/***********************************************************************//**
  118 + * @brief Construct from spatial and spectral components
  119 + *
  120 + * @param[in] spatial Spatial model component.
  121 + * @param[in] spectral Spectral model component.
  122 + *
  123 + * Constructs a background model from a @p spatial and a @p spectral model
  124 + * component. The temporal component is assumed to be constant.
  125 + ***************************************************************************/
  126 +GCTAModelBackground::GCTAModelBackground(const GCTAModelSpatial& spatial,
  127 + const GModelSpectral& spectral) :
  128 + GModelData()
  129 +{
  130 + // Initialise members
  131 + init_members();
  132 +
  133 + // Allocate temporal constant model
  134 + GModelTemporalConst temporal;
  135 +
  136 + // Clone model components
  137 + m_spatial = spatial.clone();
  138 + m_spectral = spectral.clone();
  139 + m_temporal = temporal.clone();
  140 +
  141 + // Set parameter pointers
  142 + set_pointers();
  143 +
  144 + // Return
  145 + return;
  146 +}
  147 +
  148 +
  149 +/***********************************************************************//**
  150 + * @brief Construct from model components
  151 + *
  152 + * @param[in] spatial Spatial model component.
  153 + * @param[in] spectral Spectral model component.
  154 + * @param[in] temporal Temporal model component.
  155 + *
  156 + * Constructs a background model from a @p spatial, a @p spectral and a
  157 + * @p temporal model component.
  158 + ***************************************************************************/
  159 +GCTAModelBackground::GCTAModelBackground(const GCTAModelSpatial& spatial,
  160 + const GModelSpectral& spectral,
  161 + const GModelTemporal& temporal) :
  162 + GModelData()
  163 +{
  164 + // Initialise members
  165 + init_members();
  166 +
  167 + // Clone model components
  168 + m_spatial = spatial.clone();
  169 + m_spectral = spectral.clone();
  170 + m_temporal = temporal.clone();
  171 +
  172 + // Set parameter pointers
  173 + set_pointers();
  174 +
  175 + // Return
  176 + return;
  177 +}
  178 +
  179 +
  180 +/***********************************************************************//**
  181 + * @brief Copy constructor
  182 + *
  183 + * @param[in] model Background model.
  184 + *
  185 + * Constructs a background model by copying information from an existing
  186 + * model. Note that the copy is a deep copy, so the original object can be
  187 + * destroyed after the copy without any loss of information.
  188 + ***************************************************************************/
  189 +GCTAModelBackground::GCTAModelBackground(const GCTAModelBackground& model) :
  190 + GModelData(model)
  191 +{
  192 + // Initialise private members for clean destruction
  193 + init_members();
  194 +
  195 + // Copy members
  196 + copy_members(model);
  197 +
  198 + // Set parameter pointers
  199 + set_pointers();
  200 +
  201 + // Return
  202 + return;
  203 +}
  204 +
  205 +
  206 +/***********************************************************************//**
  207 + * @brief Destructor
  208 + *
  209 + * Destroys a background model.
  210 + ***************************************************************************/
  211 +GCTAModelBackground::~GCTAModelBackground(void)
  212 +{
  213 + // Free members
  214 + free_members();
  215 +
  216 + // Return
  217 + return;
  218 +}
  219 +
  220 +
  221 +/*==========================================================================
  222 + = =
  223 + = Operators =
  224 + = =
  225 + ==========================================================================*/
  226 +
  227 +/***********************************************************************//**
  228 + * @brief Assignment operator
  229 + *
  230 + * @param[in] model Background model.
  231 + *
  232 + * Assigns the information from a background model to instance. Note that a
  233 + * deep copy of the information is performed, so the original instance can be
  234 + * destroyed after the assignment without any loss of information.
  235 + ***************************************************************************/
  236 +GCTAModelBackground& GCTAModelBackground::operator=(const GCTAModelBackground& model)
  237 +{
  238 + // Execute only if object is not identical
  239 + if (this != &model) {
  240 +
  241 + // Copy base class members
  242 + this->GModelData::operator=(model);
  243 +
  244 + // Free members
  245 + free_members();
  246 +
  247 + // Initialise members
  248 + init_members();
  249 +
  250 + // Copy members (this method also sets the parameter pointers)
  251 + copy_members(model);
  252 +
  253 + } // endif: object was not identical
  254 +
  255 + // Return
  256 + return *this;
  257 +}
  258 +
  259 +
  260 +/*==========================================================================
  261 + = =
  262 + = Public methods =
  263 + = =
  264 + ==========================================================================*/
  265 +
  266 +/***********************************************************************//**
  267 + * @brief Clear instance
  268 + *
  269 + * Resets the instance to a clean initial state. All information that resided
  270 + * in the object will be lost.
  271 + ***************************************************************************/
  272 +void GCTAModelBackground::clear(void)
  273 +{
  274 + // Free class members (base and derived classes, derived class first)
  275 + free_members();
  276 + this->GModelData::free_members();
  277 + this->GModel::free_members();
  278 +
  279 + // Initialise members
  280 + this->GModel::init_members();
  281 + this->GModelData::init_members();
  282 + init_members();
  283 +
  284 + // Return
  285 + return;
  286 +}
  287 +
  288 +
  289 +/***********************************************************************//**
  290 + * @brief Clone instance
  291 + *
  292 + * @return Pointer to deep copy of background model.
  293 + ***************************************************************************/
  294 +GCTAModelBackground* GCTAModelBackground::clone(void) const
  295 +{
  296 + return new GCTAModelBackground(*this);
  297 +}
  298 +
  299 +
  300 +/***********************************************************************//**
  301 + * @brief Set spatial model component
  302 + *
  303 + * @param[in] spatial Pointer to spatial model component.
  304 + *
  305 + * Sets the spatial model component of the model.
  306 + ***************************************************************************/
  307 +void GCTAModelBackground::spatial(const GCTAModelSpatial* spatial)
  308 +{
  309 + // Free spatial model component only if it differs from current
  310 + // component. This prevents unintentional deallocation of the argument
  311 + if ((m_spatial != NULL) && (m_spatial != spatial)) {
  312 + delete m_spatial;
  313 + }
  314 +
  315 + // Clone spatial model component if it exists, otherwise set pointer
  316 + // to NULL
  317 + m_spatial = (spatial != NULL) ? spatial->clone() : NULL;
  318 +
  319 + // Set parameter pointers
  320 + set_pointers();
  321 +
  322 + // Return
  323 + return;
  324 +}
  325 +
  326 +
  327 +/***********************************************************************//**
  328 + * @brief Set spectral model component
  329 + *
  330 + * @param[in] spectral Pointer to spectral model component.
  331 + *
  332 + * Sets the spectral model component of the model.
  333 + ***************************************************************************/
  334 +void GCTAModelBackground::spectral(const GModelSpectral* spectral)
  335 +{
  336 + // Free spectral model component only if it differs from current
  337 + // component. This prevents unintential deallocation of the argument
  338 + if ((m_spectral != NULL) && (m_spectral != spectral)) {
  339 + delete m_spectral;
  340 + }
  341 +
  342 + // Clone spectral model component if it exists, otherwise set pointer
  343 + // to NULL
  344 + m_spectral = (spectral != NULL) ? spectral->clone() : NULL;
  345 +
  346 + // Set parameter pointers
  347 + set_pointers();
  348 +
  349 + // Return
  350 + return;
  351 +}
  352 +
  353 +
  354 +/***********************************************************************//**
  355 + * @brief Set temporal model component
  356 + *
  357 + * @param[in] temporal Pointer to temporal model component.
  358 + *
  359 + * Sets the temporal model component of the model.
  360 + ***************************************************************************/
  361 +void GCTAModelBackground::temporal(const GModelTemporal* temporal)
  362 +{
  363 + // Free temporal model component only if it differs from current
  364 + // component. This prevents unintential deallocation of the argument
  365 + if ((m_temporal != NULL) && (m_temporal != temporal)) {
  366 + delete m_temporal;
  367 + }
  368 +
  369 + // Clone temporal model component if it exists, otherwise set pointer
  370 + // to NULL
  371 + m_temporal = (temporal != NULL) ? temporal->clone() : NULL;
  372 +
  373 + // Set parameter pointers
  374 + set_pointers();
  375 +
  376 + // Return
  377 + return;
  378 +}
  379 +
  380 +
  381 +/***********************************************************************//**
  382 + * @brief Evaluate function
  383 + *
  384 + * @param[in] event Observed event.
  385 + * @param[in] obs Observation.
  386 + * @param[in] gradients Compute gradients?
  387 + * @return Function value.
  388 + *
  389 + * Evaluates the background model which is a factorization of a spatial,
  390 + * spectral and temporal model component. This method also applies a deadtime
  391 + * correction factor, so that the normalization of the model is a real rate
  392 + * (counts/exposure time).
  393 + *
  394 + * If the @p gradients flag is true the method will also set the parameter
  395 + * gradients of the model parameters.
  396 + *
  397 + * @todo Add bookkeeping of last value and evaluate only if argument
  398 + * changed
  399 + * @todo Verify that CTA instrument direction pointer is valid, or better,
  400 + * add an offset method to GCTAPointing. Ideally, we should precompute
  401 + * all offset angles (for an event cube this may only be done easily
  402 + * if the pointing has been fixed; otherwise we need a structure
  403 + * similar to the Fermi/LAT livetime cube that provides the effective
  404 + * sky exposure as function of offset angle).
  405 + ***************************************************************************/
  406 +double GCTAModelBackground::eval(const GEvent& event,
  407 + const GObservation& obs,
  408 + const bool& gradients) const
  409 +{
  410 + // Get reference on CTA pointing from observation and reference on CTA
  411 + // instrument direction from event
  412 + const GCTAPointing& pnt = gammalib::cta_pnt(G_EVAL, obs);
  413 + const GCTAInstDir& dir = gammalib::cta_dir(G_EVAL, event);
  414 +
  415 + // Compute offset angle (in degrees)
  416 + //double offset = dir.dir().dist_deg(pnt.dir());
  417 +
  418 + // Evaluate function and gradients
  419 + double spat = (spatial() != NULL)
  420 + ? spatial()->eval(dir, event.energy(), event.time(), gradients) : 1.0;
  421 + double spec = (spectral() != NULL)
  422 + ? spectral()->eval(event.energy(), event.time(), gradients)
  423 + : 1.0;
  424 + double temp = (temporal() != NULL)
  425 + ? temporal()->eval(event.time(), gradients) : 1.0;
  426 +
  427 + // Compute value
  428 + double value = spat * spec * temp;
  429 +
  430 + // Apply deadtime correction
  431 + double deadc = obs.deadc(event.time());
  432 + value *= deadc;
  433 +
  434 + // Optionally compute partial derivatives
  435 + if (gradients) {
  436 +
  437 + // Multiply factors to radial gradients
  438 + if (spatial() != NULL) {
  439 + double fact = spec * temp * deadc;
  440 + if (fact != 1.0) {
  441 + for (int i = 0; i < spatial()->size(); ++i)
  442 + (*spatial())[i].factor_gradient((*spatial())[i].factor_gradient() * fact);
  443 + }
  444 + }
  445 +
  446 + // Multiply factors to spectral gradients
  447 + if (spectral() != NULL) {
  448 + double fact = spat * temp * deadc;
  449 + if (fact != 1.0) {
  450 + for (int i = 0; i < spectral()->size(); ++i)
  451 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
  452 + }
  453 + }
  454 +
  455 + // Multiply factors to temporal gradients
  456 + if (temporal() != NULL) {
  457 + double fact = spat * spec * deadc;
  458 + if (fact != 1.0) {
  459 + for (int i = 0; i < temporal()->size(); ++i)
  460 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
  461 + }
  462 + }
  463 +
  464 + } // endif: computed partial derivatives
  465 +
  466 + // Return
  467 + return value;
  468 +}
  469 +
  470 +
  471 +/***********************************************************************//**
  472 + * @brief Return spatially integrated data model
  473 + *
  474 + * @param[in] obsEng Measured event energy.
  475 + * @param[in] obsTime Measured event time.
  476 + * @param[in] obs Observation.
  477 + *
  478 + * @exception GException::no_list
  479 + * No valid CTA event list found in observation
  480 + * @exception GCTAException::no_pointing
  481 + * No valid CTA pointing found in observation
  482 + *
  483 + * Spatially integrates the data model for a given measured event energy and
  484 + * event time. This method also applies a deadtime correction factor, so that
  485 + * the normalization of the model is a real rate (counts/exposure time).
  486 + ***************************************************************************/
  487 +double GCTAModelBackground::npred(const GEnergy& obsEng,
  488 + const GTime& obsTime,
  489 + const GObservation& obs) const
  490 +{
  491 + // Initialise result
  492 + double npred = 0.0;
  493 +
  494 + // Evaluate only if model is valid
  495 + if (valid_model()) {
  496 +
  497 + // Retrieve CTA pointing and event list
  498 + const GCTAPointing& pnt = gammalib::cta_pnt(G_NPRED, obs);
  499 + const GCTAEventList& events = gammalib::cta_event_list(G_NPRED, obs);
  500 +
  501 + // Get ROI radius in radians
  502 + double roi_radius = events.roi().radius() * gammalib::deg2rad;
  503 +
  504 + // Get distance from ROI centre in radians
  505 + double roi_distance = events.roi().centre().dir().dist(pnt.dir());
  506 +
  507 + // Setup integration function
  508 + /*
  509 + GCTAModelBackground::roi_kern integrand(radial(), roi_radius, roi_distance);
  510 +
  511 + // Setup integrator
  512 + GIntegral integral(&integrand);
  513 +
  514 + // Setup integration boundaries
  515 + double rmin = (roi_distance > roi_radius) ? roi_distance-roi_radius : 0.0;
  516 + double rmax = roi_radius + roi_distance;
  517 +
  518 + // Spatially integrate radial component
  519 + npred = integral.romberg(rmin, rmax);
  520 +
  521 + // Multiply in spectral and temporal components
  522 + npred *= spectral()->eval(obsEng, obsTime);
  523 + npred *= temporal()->eval(obsTime);
  524 +
  525 + // Apply deadtime correction
  526 + npred *= obs.deadc(obsTime);
  527 + */
  528 +
  529 + } // endif: model was valid
  530 +
  531 + // Return
  532 + return npred;
  533 +}
  534 +
  535 +
  536 +/***********************************************************************//**
  537 + * @brief Return simulated list of events
  538 + *
  539 + * @param[in] obs Observation.
  540 + * @param[in] ran Random number generator.
  541 + *
  542 + * @exception GCTAException::no_pointing
  543 + * No CTA pointing found in observation.
  544 + *
  545 + * Draws a sample of events from the background model using a Monte Carlo
  546 + * simulation. The pointing information, the energy boundaries and the good
  547 + * time interval for the sampling will be extracted from the observation
  548 + * argument that is passed to the method. The method also requires a random
  549 + * number generator of type GRan which is passed by reference, hence the
  550 + * state of the random number generator will be changed by the method.
  551 + *
  552 + * The method also applies a deadtime correction using a Monte Carlo process,
  553 + * taking into account temporal deadtime variations. For this purpose, the
  554 + * method makes use of the time dependent GObservation::deadc method.
  555 + ***************************************************************************/
  556 +GCTAEventList* GCTAModelBackground::mc(const GObservation& obs,
  557 + GRan& ran) const
  558 +{
  559 + // Initialise new event list
  560 + GCTAEventList* list = new GCTAEventList;
  561 +
  562 + // Continue only if model is valid)
  563 + if (valid_model()) {
  564 +
  565 + // Retrieve CTA pointing
  566 + const GCTAPointing& pnt = gammalib::cta_pnt(G_NPRED, obs);
  567 +
  568 + // Convert CTA pointing direction in instrument system
  569 + GCTAInstDir pnt_dir(pnt.dir());
  570 +
  571 + // Loop over all energy boundaries
  572 + for (int ieng = 0; ieng < obs.events()->ebounds().size(); ++ieng) {
  573 +
  574 + // Compute the on-axis background rate in model within the
  575 + // energy boundaries from spectral component (units: cts/s/sr)
  576 + double flux = spectral()->flux(obs.events()->ebounds().emin(ieng),
  577 + obs.events()->ebounds().emax(ieng));
  578 +
  579 + // Compute solid angle used for normalization
  580 + double area = spatial()->omega(); //TODO: Check if this is reasonable
  581 +
  582 + // Derive expecting rate (units: cts/s). Note that the time here
  583 + // is good time. Deadtime correction will be done later.
  584 + double rate = flux * area;
  585 +
  586 + // Debug option: dump rate
  587 + #if defined(G_DUMP_MC)
  588 + std::cout << "GCTAModelBackground::mc(\"" << name() << "\": ";
  589 + std::cout << "flux=" << flux << " cts/s/sr, ";
  590 + std::cout << "area=" << area << " sr, ";
  591 + std::cout << "rate=" << rate << " cts/s)" << std::endl;
  592 + #endif
  593 +
  594 + // Loop over all good time intervals
  595 + for (int itime = 0; itime < obs.events()->gti().size(); ++itime) {
  596 +
  597 + // Get event arrival times from temporal model
  598 + GTimes times = m_temporal->mc(rate,
  599 + obs.events()->gti().tstart(itime),
  600 + obs.events()->gti().tstop(itime),
  601 + ran);
  602 +
  603 + // Get number of events
  604 + int n_events = times.size();
  605 +
  606 + // Reserve space for events
  607 + if (n_events > 0) {
  608 + list->reserve(n_events);
  609 + }
  610 +
  611 + // Loop over events
  612 + for (int i = 0; i < n_events; ++i) {
  613 +
  614 + // Apply deadtime correction
  615 + double deadc = obs.deadc(times[i]);
  616 + if (deadc < 1.0) {
  617 + if (ran.uniform() > deadc) {
  618 + continue;
  619 + }
  620 + }
  621 +
  622 + // Set event energy
  623 + GEnergy energy = spectral()->mc(obs.events()->ebounds().emin(ieng),
  624 + obs.events()->ebounds().emax(ieng),
  625 + times[i],
  626 + ran);
  627 +
  628 + // Set event direction (only DETX and DETY are set)
  629 + GCTAInstDir dir = spatial()->mc(energy, times[i], ran);
  630 +
  631 + // Compute RA/DEC
  632 + //TODO
  633 +
  634 + // Allocate event
  635 + GCTAEventAtom event;
  636 +
  637 + // Set event attributes
  638 + event.dir(dir);
  639 + event.energy(energy);
  640 + event.time(times[i]);
  641 +
  642 + // Append event to list
  643 + list->append(event);
  644 +
  645 + } // endfor: looped over all events
  646 +
  647 + } // endfor: looped over all GTIs
  648 +
  649 + } // endfor: looped over all energy boundaries
  650 +
  651 + } // endif: model was valid
  652 +
  653 + // Return
  654 + return list;
  655 +}
  656 +
  657 +
  658 +/***********************************************************************//**
  659 + * @brief Read model from XML element
  660 + *
  661 + * @param[in] xml XML element.
  662 + *
  663 + * Reads the sky model from an XML element. The XML element is expected to
  664 + * respect the following format:
  665 + *
  666 + * <source name=".." type=".." instrument=".." id="..">
  667 + * <spectrum type="..">
  668 + * ..
  669 + * </spectrum>
  670 + * <spatialModel type="..">
  671 + * ..
  672 + * </spatialModel>
  673 + * <temporal type="..">
  674 + * ..
  675 + * </temporal>
  676 + * </source>
  677 + *
  678 + * The temporal element is optional. In no temporal element is specified a
  679 + * constant component with unity normalization will be assumed.
  680 + ***************************************************************************/
  681 +void GCTAModelBackground::read(const GXmlElement& xml)
  682 +{
  683 + // Clear model
  684 + clear();
  685 +
  686 + // Get pointers on spatial and spectral model components
  687 + const GXmlElement* spat = xml.element("spatialModel", 0);
  688 + const GXmlElement* spec = xml.element("spectrum", 0);
  689 +
  690 + // Set spatial and spectral model components
  691 + m_spatial = xml_spatial(*spat);
  692 + m_spectral = xml_spectral(*spec);
  693 +
  694 + // Handle optional temporal model
  695 + if (xml.elements("temporal") > 0) {
  696 + const GXmlElement* temp = xml.element("temporal", 0);
  697 + m_temporal = xml_temporal(*temp);
  698 + }
  699 + else {
  700 + GModelTemporalConst temporal;
  701 + m_temporal = temporal.clone();
  702 + }
  703 +
  704 + // Read model attributes
  705 + read_attributes(xml);
  706 +
  707 + // Set parameter pointers
  708 + set_pointers();
  709 +
  710 + // Return
  711 + return;
  712 +}
  713 +
  714 +
  715 +/***********************************************************************//**
  716 + * @brief Write model into XML element
  717 + *
  718 + * @param[in] xml XML element.
  719 + *
  720 + * Writes the sky model into an XML source library. The format written to
  721 + * the @p xml element is as follows:
  722 + *
  723 + * <source name=".." type=".." instrument=".." id="..">
  724 + * <spectrum type="..">
  725 + * ..
  726 + * </spectrum>
  727 + * <spatialModel type="..">
  728 + * ..
  729 + * </spatialModel>
  730 + * <temporal type="..">
  731 + * ..
  732 + * </temporal>
  733 + * </source>
  734 + *
  735 + * For compatibility reasons the temporal element will only be written if it
  736 + * is a non-constant component or a constant component with a normalization
  737 + * that differs from unity.
  738 + ***************************************************************************/
  739 +void GCTAModelBackground::write(GXmlElement& xml) const
  740 +{
  741 + // Initialise pointer on source
  742 + GXmlElement* src = NULL;
  743 +
  744 + // Search corresponding source
  745 + int n = xml.elements("source");
  746 + for (int k = 0; k < n; ++k) {
  747 + GXmlElement* element = xml.element("source", k);
  748 + if (element->attribute("name") == name()) {
  749 + src = element;
  750 + break;
  751 + }
  752 + }
  753 +
  754 + // If the temporal model is not a constant with unit normalization then
  755 + // set cons to a NULL pointer
  756 + GModelTemporalConst* cons = dynamic_cast<GModelTemporalConst*>(temporal());
  757 + if (cons != NULL) {
  758 + if (cons->norm() != 1.0) {
  759 + cons = NULL;
  760 + }
  761 + }
  762 +
  763 + // If no source with corresponding name was found then append one
  764 + if (src == NULL) {
  765 + src = xml.append("source");
  766 + if (spatial() != NULL) src->append(GXmlElement("spatialModel"));
  767 + if (spectral() != NULL) src->append(GXmlElement("spectrum"));
  768 + if (temporal() != NULL && cons == NULL) src->append(GXmlElement("temporal"));
  769 + }
  770 +
  771 + // Write spatial model
  772 + if (spatial() != NULL) {
  773 + GXmlElement* spat = src->element("spatialModel", 0);
  774 + spatial()->write(*spat);
  775 + }
  776 +
  777 + // Write spectral model
  778 + if (spectral() != NULL) {
  779 + GXmlElement* spec = src->element("spectrum", 0);
  780 + spectral()->write(*spec);
  781 + }
  782 +
  783 + // Write temporal model (only if not a constant with unit normalization
  784 + // factor)
  785 + if (temporal() != NULL && cons == NULL) {
  786 + GXmlElement* temp = src->element("temporal", 0);
  787 + temporal()->write(*temp);
  788 + }
  789 +
  790 + // Write model attributes
  791 + write_attributes(*src);
  792 +
  793 + // Return
  794 + return;
  795 +}
  796 +
  797 +
  798 +/***********************************************************************//**
  799 + * @brief Print model information
  800 + *
  801 + * @param[in] chatter Chattiness.
  802 + * @return String containing model information.
  803 + ***************************************************************************/
  804 +std::string GCTAModelBackground::print(const GChatter& chatter) const
  805 +{
  806 + // Initialise result string
  807 + std::string result;
  808 +
  809 + // Continue only if chatter is not silent
  810 + if (chatter != SILENT) {
  811 +
  812 + // Append header
  813 + result.append("=== GCTAModelBackground ===");
  814 +
  815 + // Determine number of parameters per type
  816 + int n_spatial = (spatial() != NULL) ? spatial()->size() : 0;
  817 + int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
  818 + int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
  819 +
  820 + // Append attributes
  821 + result.append("\n"+print_attributes());
  822 +
  823 + // Append model type
  824 + result.append("\n"+gammalib::parformat("Model type")+type());
  825 +
  826 + // Append model type
  827 + result.append("\n"+gammalib::parformat("Model components"));
  828 + if (n_spatial > 0) {
  829 + result.append("\""+spatial()->type()+"\"");
  830 + if (n_spectral > 0 || n_temporal > 0) {
  831 + result.append(" * ");
  832 + }
  833 + }
  834 + if (n_spectral > 0) {
  835 + result.append("\""+spectral()->type()+"\"");
  836 + if (n_temporal > 0) {
  837 + result.append(" * ");
  838 + }
  839 + }
  840 + if (n_temporal > 0) {
  841 + result.append("\""+temporal()->type()+"\"");
  842 + }
  843 +
  844 + // Append parameters
  845 + result.append("\n"+gammalib::parformat("Number of parameters") +
  846 + gammalib::str(size()));
  847 + result.append("\n"+gammalib::parformat("Number of spatial par's") +
  848 + gammalib::str(n_spatial));
  849 + for (int i = 0; i < n_spatial; ++i) {
  850 + result.append("\n"+(*spatial())[i].print());
  851 + }
  852 + result.append("\n"+gammalib::parformat("Number of spectral par's") +
  853 + gammalib::str(n_spectral));
  854 + for (int i = 0; i < n_spectral; ++i) {
  855 + result.append("\n"+(*spectral())[i].print());
  856 + }
  857 + result.append("\n"+gammalib::parformat("Number of temporal par's") +
  858 + gammalib::str(n_temporal));
  859 + for (int i = 0; i < n_temporal; ++i) {
  860 + result.append("\n"+(*temporal())[i].print());
  861 + }
  862 +
  863 + } // endif: chatter was not silent
  864 +
  865 + // Return result
  866 + return result;
  867 +}
  868 +
  869 +
  870 +/*==========================================================================
  871 + = =
  872 + = Private methods =
  873 + = =
  874 + ==========================================================================*/
  875 +
  876 +/***********************************************************************//**
  877 + * @brief Initialise class members
  878 + ***************************************************************************/
  879 +void GCTAModelBackground::init_members(void)
  880 +{
  881 + // Initialise members
  882 + m_spatial = NULL;
  883 + m_spectral = NULL;
  884 + m_temporal = NULL;
  885 +
  886 + // Return
  887 + return;
  888 +}
  889 +
  890 +
  891 +/***********************************************************************//**
  892 + * @brief Copy class members
  893 + *
  894 + * @param[in] model Background model.
  895 + ***************************************************************************/
  896 +void GCTAModelBackground::copy_members(const GCTAModelBackground& model)
  897 +{
  898 + // Clone radial, spectral and temporal model components
  899 + m_spatial = (model.m_spatial != NULL) ? model.m_spatial->clone() : NULL;
  900 + m_spectral = (model.m_spectral != NULL) ? model.m_spectral->clone() : NULL;
  901 + m_temporal = (model.m_temporal != NULL) ? model.m_temporal->clone() : NULL;
  902 +
  903 + // Set parameter pointers
  904 + set_pointers();
  905 +
  906 + // Return
  907 + return;
  908 +}
  909 +
  910 +
  911 +/***********************************************************************//**
  912 + * @brief Delete class members
  913 + ***************************************************************************/
  914 +void GCTAModelBackground::free_members(void)
  915 +{
  916 + // Free memory
  917 + if (m_spatial != NULL) delete m_spatial;
  918 + if (m_spectral != NULL) delete m_spectral;
  919 + if (m_temporal != NULL) delete m_temporal;
  920 +
  921 + // Signal free pointers
  922 + m_spatial = NULL;
  923 + m_spectral = NULL;
  924 + m_temporal = NULL;
  925 +
  926 + // Return
  927 + return;
  928 +}
  929 +
  930 +
  931 +/***********************************************************************//**
  932 + * @brief Set pointers
  933 + *
  934 + * Set pointers to all model parameters. The pointers are stored in a vector
  935 + * that is member of the GModelData base class.
  936 + ***************************************************************************/
  937 +void GCTAModelBackground::set_pointers(void)
  938 +{
  939 + // Clear parameters
  940 + m_pars.clear();
  941 +
  942 + // Determine the number of parameters
  943 + int n_spatial = (spatial() != NULL) ? spatial()->size() : 0;
  944 + int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
  945 + int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
  946 + int n_pars = n_spatial + n_spectral + n_temporal;
  947 +
  948 + // Continue only if there are parameters
  949 + if (n_pars > 0) {
  950 +
  951 + // Gather spatial parameter pointers
  952 + for (int i = 0; i < n_spatial; ++i) {
  953 + m_pars.push_back(&((*spatial())[i]));
  954 + }
  955 +
  956 + // Gather spectral parameters
  957 + for (int i = 0; i < n_spectral; ++i) {
  958 + m_pars.push_back(&((*spectral())[i]));
  959 + }
  960 +
  961 + // Gather temporal parameters
  962 + for (int i = 0; i < n_temporal; ++i) {
  963 + m_pars.push_back(&((*temporal())[i]));
  964 + }
  965 +
  966 + }
  967 +
  968 + // Return
  969 + return;
  970 +}
  971 +
  972 +
  973 +/***********************************************************************//**
  974 + * @brief Return pointer to spatial model from XML element
  975 + *
  976 + * @param[in] spatial XML element.
  977 + * @return Pointer to spatial model.
  978 + *
  979 + * Returns pointer to a radial model that is defined in an XML element.
  980 + ***************************************************************************/
  981 +GCTAModelSpatial* GCTAModelBackground::xml_spatial(const GXmlElement& spatial) const
  982 +{
  983 + // Get radial model
  984 + GCTAModelSpatialRegistry registry;
  985 + GCTAModelSpatial* ptr = registry.alloc(spatial);
  986 +
  987 + // Return pointer
  988 + return ptr;
  989 +}
  990 +
  991 +
  992 +/***********************************************************************//**
  993 + * @brief Return pointer to spectral model from XML element
  994 + *
  995 + * @param[in] spectral XML element.
  996 + * @return Pointer to spectral model.
  997 + *
  998 + * Returns pointer to spectral model that is defined in an XML element.
  999 + ***************************************************************************/
  1000 +GModelSpectral* GCTAModelBackground::xml_spectral(const GXmlElement& spectral) const
  1001 +{
  1002 + // Get spectral model
  1003 + GModelSpectralRegistry registry;
  1004 + GModelSpectral* ptr = registry.alloc(spectral);
  1005 +
  1006 + // Return pointer
  1007 + return ptr;
  1008 +}
  1009 +
  1010 +
  1011 +/***********************************************************************//**
  1012 + * @brief Return pointer to temporal model from XML element
  1013 + *
  1014 + * @param[in] temporal XML element.
  1015 + * @return Pointer to temporal model.
  1016 + *
  1017 + * Returns pointer to temporal model that is defined in an XML element.
  1018 + ***************************************************************************/
  1019 +GModelTemporal* GCTAModelBackground::xml_temporal(const GXmlElement& temporal) const
  1020 +{
  1021 + // Get temporal model
  1022 + GModelTemporalRegistry registry;
  1023 + GModelTemporal* ptr = registry.alloc(temporal);
  1024 +
  1025 + // Return pointer
  1026 + return ptr;
  1027 +}
  1028 +
  1029 +
  1030 +/***********************************************************************//**
  1031 + * @brief Verifies if model has all components
  1032 + *
  1033 + * Returns 'true' if models has a spatial, a spectral and a temporal
  1034 + * component. Otherwise returns 'false'.
  1035 + ***************************************************************************/
  1036 +bool GCTAModelBackground::valid_model(void) const
  1037 +{
  1038 + // Set result
  1039 + bool result = ((spatial() != NULL) &&
  1040 + (spectral() != NULL) &&
  1041 + (temporal() != NULL));
  1042 +
  1043 + // Return result
  1044 + return result;
  1045 +}
  1046 +
  1047 +
  1048 +/***********************************************************************//**
  1049 + * @brief Integration kernel for the Npred method
  1050 + *
  1051 + * @param[in] offset Offset angle (radians).
  1052 + *
  1053 + * Computes
  1054 + * \f[N_{\rm pred} = \int_{\rm ROI} r(\theta) \sin \theta {\rm d}\phi\f]
  1055 + * where
  1056 + * \f$r(\theta)\f$ is the radial acceptance function,
  1057 + * \f$\theta\f$ is the measured offset angle from the pointing direction
  1058 + * in radians, and
  1059 + * \f$\phi\f$ is the measured azimuth angle. The integration is done over
  1060 + * the arc of the azimuth angle that lies within the ROI. This integration
  1061 + * is done analytically using the "cta_roi_arclength" support function.
  1062 + ***************************************************************************/
  1063 +double GCTAModelBackground::roi_kern::eval(const double& offset)
  1064 +{
  1065 + // Circumvent const correctness
  1066 + GCTAModelSpatial* spatial = const_cast<GCTAModelSpatial*>(m_parent);
  1067 +
  1068 + // Initialise Npred value
  1069 + double value = 0.0;
  1070 +
  1071 + /*
  1072 + // Continue only if offset > 0
  1073 + if (offset > 0.0) {
  1074 +
  1075 + // Get arclength for given radius in radians.
  1076 + double phi = gammalib::cta_roi_arclength(offset,
  1077 + m_dist,
  1078 + m_cosdist,
  1079 + m_sindist,
  1080 + m_roi,
  1081 + m_cosroi);
  1082 +
  1083 + // Get kernel value if phi > 0
  1084 + if (phi > 0.0) {
  1085 + value = radial->eval(offset*gammalib::rad2deg) * phi *
  1086 + std::sin(offset);
  1087 + }
  1088 +
  1089 + } // endif: offset was positive
  1090 + */
  1091 +
  1092 + // Return
  1093 + return value;
  1094 +}
... ...
inst/cta/src/GCTAModelRadial.cpp
  View file @513f99b
1 1 /***************************************************************************
2 2 * GCTAModelRadial.cpp - Abstract radial model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011 by Jurgen Knodlseder *
  4 + * copyright (C) 2011-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -21,7 +21,7 @@
21 21 /**
22 22 * @file GCTAModelRadial.cpp
23 23 * @brief Abstract radial acceptance model class implementation
24   - * @author J. Knodlseder
  24 + * @author Juergen Knoedlseder
25 25 */
26 26  
27 27 /* __ Includes ___________________________________________________________ */
... ... @@ -31,6 +31,7 @@
31 31 #include "GException.hpp"
32 32 #include "GTools.hpp"
33 33 #include "GCTAModelRadial.hpp"
  34 +#include "GCTAInstDir.hpp"
34 35  
35 36 /* __ Method name definitions ____________________________________________ */
36 37 #define G_ACCESS1 "GModelSpatial::operator[](int&)"
... ... @@ -225,6 +226,46 @@ const GModelPar&amp; GCTAModelRadial::operator[](const std::string&amp; name) const
225 226 = =
226 227 ==========================================================================*/
227 228  
  229 +/***********************************************************************//**
  230 + * @brief Evaluate function
  231 + *
  232 + * @param[in] dir Event direction.
  233 + * @param[in] energy Event energy.
  234 + * @param[in] time Event time.
  235 + * @param[in] gradients Compute gradients?
  236 + * @return Function value
  237 + *
  238 + * @todo Implement method
  239 + ***************************************************************************/
  240 +double GCTAModelRadial::eval(const GCTAInstDir& dir,
  241 + const GEnergy& energy,
  242 + const GTime& time,
  243 + const bool& gradients) const
  244 +{
  245 + // Return
  246 + return 0.0;
  247 +}
  248 +
  249 +
  250 +/***********************************************************************//**
  251 + * @brief Returns MC instrument direction
  252 + *
  253 + * @param[in] energy Event energy.
  254 + * @param[in] time Event time.
  255 + * @param[in,out] ran Random number generator.
  256 + * @return Instrument direction
  257 + *
  258 + * @todo Implement method
  259 + ***************************************************************************/
  260 +GCTAInstDir GCTAModelRadial::mc(const GEnergy& energy,
  261 + const GTime& time,
  262 + GRan& ran) const
  263 +{
  264 + // Return
  265 + return GCTAInstDir();
  266 +}
  267 +
  268 +
228 269 /*==========================================================================
229 270 = =
230 271 = Private methods =
... ...
inst/cta/src/GCTAModelRadialGauss.cpp
  View file @513f99b
1 1 /***************************************************************************
2 2 * GCTAModelRadialGauss.cpp - Radial Gaussian CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -32,15 +32,18 @@
32 32 #include "GException.hpp"
33 33 #include "GTools.hpp"
34 34 #include "GMath.hpp"
  35 +#include "GRan.hpp"
35 36 #include "GIntegral.hpp"
36 37 #include "GCTAModelRadialGauss.hpp"
37 38 #include "GCTAModelRadialRegistry.hpp"
  39 +#include "GCTAModelSpatialRegistry.hpp"
38 40  
39 41 /* __ Constants __________________________________________________________ */
40 42  
41 43 /* __ Globals ____________________________________________________________ */
42   -const GCTAModelRadialGauss g_cta_radial_gauss_seed;
43   -const GCTAModelRadialRegistry g_cta_radial_gauss_registry(&g_cta_radial_gauss_seed);
  44 +const GCTAModelRadialGauss g_cta_radial_gauss_seed;
  45 +const GCTAModelRadialRegistry g_cta_radial_gauss_registry(&g_cta_radial_gauss_seed);
  46 +const GCTAModelSpatialRegistry g_cta_radial_gauss_spatial_registry(&g_cta_radial_gauss_seed);
44 47  
45 48 /* __ Method name definitions ____________________________________________ */
46 49 #define G_READ "GCTAModelRadialGauss::read(GXmlElement&)"
... ...
inst/cta/src/GCTAModelRadialPolynom.cpp
  View file @513f99b
1 1 /***************************************************************************
2 2 * GCTAModelRadialPolynom.cpp - Radial Polynom CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -32,16 +32,19 @@
32 32 #include "GException.hpp"
33 33 #include "GTools.hpp"
34 34 #include "GMath.hpp"
  35 +#include "GRan.hpp"
35 36 #include "GIntegral.hpp"
36 37 #include "GCTAModelRadialPolynom.hpp"
37 38 #include "GCTAModelRadialRegistry.hpp"
  39 +#include "GCTAModelSpatialRegistry.hpp"
38 40  
39 41 /* __ Constants __________________________________________________________ */
40 42 const double g_cta_radial_polynom_offset_max = 4.0;
41 43  
42 44 /* __ Globals ____________________________________________________________ */
43   -const GCTAModelRadialPolynom g_cta_radial_polynom_seed;
44   -const GCTAModelRadialRegistry g_cta_radial_polynom_registry(&g_cta_radial_polynom_seed);
  45 +const GCTAModelRadialPolynom g_cta_radial_polynom_seed;
  46 +const GCTAModelRadialRegistry g_cta_radial_polynom_registry(&g_cta_radial_polynom_seed);
  47 +const GCTAModelSpatialRegistry g_cta_radial_polynom_spatial_registry(&g_cta_radial_polynom_seed);
45 48  
46 49 /* __ Method name definitions ____________________________________________ */
47 50 #define G_READ "GCTAModelRadialPolynom::read(GXmlElement&)"
... ...
inst/cta/src/GCTAModelRadialProfile.cpp
  View file @513f99b
1 1 /***************************************************************************
2 2 * GCTAModelRadialProfile.cpp - Radial Profile CTA model class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -32,15 +32,18 @@
32 32 #include "GException.hpp"
33 33 #include "GTools.hpp"
34 34 #include "GMath.hpp"
  35 +#include "GRan.hpp"
35 36 #include "GIntegral.hpp"
36 37 #include "GCTAModelRadialProfile.hpp"
37 38 #include "GCTAModelRadialRegistry.hpp"
  39 +#include "GCTAModelSpatialRegistry.hpp"
38 40  
39 41 /* __ Constants __________________________________________________________ */
40 42  
41 43 /* __ Globals ____________________________________________________________ */
42   -const GCTAModelRadialProfile g_cta_radial_profile_seed;
43   -const GCTAModelRadialRegistry g_cta_radial_profile_registry(&g_cta_radial_profile_seed);
  44 +const GCTAModelRadialProfile g_cta_radial_profile_seed;
  45 +const GCTAModelRadialRegistry g_cta_radial_profile_registry(&g_cta_radial_profile_seed);
  46 +const GCTAModelSpatialRegistry g_cta_radial_profile_spatial_registry(&g_cta_radial_profile_seed);
44 47  
45 48 /* __ Method name definitions ____________________________________________ */
46 49 #define G_READ "GCTAModelRadialProfile::read(GXmlElement&)"
... ...
inst/cta/src/GCTAModelSpatial.cpp 0 → 100644
  View file @513f99b
  1 +/***************************************************************************
  2 + * GCTAModelSpatial.cpp - Spatial model abstract base class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Jurgen Knodlseder *
  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 GCTAModelSpatial.cpp
  23 + * @brief Abstract spatial model class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GTools.hpp"
  33 +#include "GCTAModelSpatial.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_ACCESS1 "GModelSpatial::operator[](int&)"
  37 +#define G_ACCESS2 "GModelSpatial::operator[](std::string&)"
  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 +GCTAModelSpatial::GCTAModelSpatial(void)
  56 +{
  57 + // Initialise members
  58 + init_members();
  59 +
  60 + // Return
  61 + return;
  62 +}
  63 +
  64 +
  65 +/***********************************************************************//**
  66 + * @brief Copy constructor
  67 + *
  68 + * @param[in] model Spatial model.
  69 + ***************************************************************************/
  70 +GCTAModelSpatial::GCTAModelSpatial(const GCTAModelSpatial& model)
  71 +{
  72 + // Initialise members
  73 + init_members();
  74 +
  75 + // Copy members
  76 + copy_members(model);
  77 +
  78 + // Return
  79 + return;
  80 +}
  81 +
  82 +
  83 +/***********************************************************************//**
  84 + * @brief Destructor
  85 + ***************************************************************************/
  86 +GCTAModelSpatial::~GCTAModelSpatial(void)
  87 +{
  88 + // Free members
  89 + free_members();
  90 +
  91 + // Return
  92 + return;
  93 +}
  94 +
  95 +
  96 +/*==========================================================================
  97 + = =
  98 + = Operators =
  99 + = =
  100 + ==========================================================================*/
  101 +
  102 +/***********************************************************************//**
  103 + * @brief Assignment operator
  104 + *
  105 + * @param[in] model Spatial model.
  106 + ***************************************************************************/
  107 +GCTAModelSpatial& GCTAModelSpatial::operator=(const GCTAModelSpatial& model)
  108 +{
  109 + // Execute only if object is not identical
  110 + if (this != &model) {
  111 +
  112 + // Free members
  113 + free_members();
  114 +
  115 + // Initialise members
  116 + init_members();
  117 +
  118 + // Copy members
  119 + copy_members(model);
  120 +
  121 + } // endif: object was not identical
  122 +
  123 + // Return
  124 + return *this;
  125 +}
  126 +
  127 +
  128 +/***********************************************************************//**
  129 + * @brief Returns model parameter
  130 + *
  131 + * @param[in] index Parameter index [0,...,size()-1].
  132 + *
  133 + * @exception GException::out_of_range
  134 + * Parameter index is out of range.
  135 + ***************************************************************************/
  136 +GModelPar& GCTAModelSpatial::operator[](const int& index)
  137 +{
  138 + // Compile option: raise exception if index is out of range
  139 + #if defined(G_RANGE_CHECK)
  140 + if (index < 0 || index >= size()) {
  141 + throw GException::out_of_range(G_ACCESS1, "Spatial parameter index",
  142 + index, size());
  143 + }
  144 + #endif
  145 +
  146 + // Return reference
  147 + return *(m_pars[index]);
  148 +}
  149 +
  150 +
  151 +/***********************************************************************//**
  152 + * @brief Returns model parameter (const version)
  153 + *
  154 + * @param[in] index Parameter index [0,...,size()-1].
  155 + *
  156 + * @exception GException::out_of_range
  157 + * Parameter index is out of range.
  158 + ***************************************************************************/
  159 +const GModelPar& GCTAModelSpatial::operator[](const int& index) const
  160 +{
  161 + // Compile option: raise exception if index is out of range
  162 + #if defined(G_RANGE_CHECK)
  163 + if (index < 0 || index >= size()) {
  164 + throw GException::out_of_range(G_ACCESS1, "Spatial parameter index",
  165 + index, size());
  166 + }
  167 + #endif
  168 +
  169 + // Return reference
  170 + return *(m_pars[index]);
  171 +}
  172 +
  173 +
  174 +/***********************************************************************//**
  175 + * @brief Returns reference to model parameter
  176 + *
  177 + * @param[in] name Parameter name.
  178 + *
  179 + * @exception GException::invalid_argument
  180 + * Parameter with specified name not found in container.
  181 + ***************************************************************************/
  182 +GModelPar& GCTAModelSpatial::operator[](const std::string& name)
  183 +{
  184 + // Get parameter index
  185 + int index = 0;
  186 + for (; index < size(); ++index) {
  187 + if (m_pars[index]->name() == name)
  188 + break;
  189 + }
  190 +
  191 + // Throw exception if parameter name was not found
  192 + if (index >= size()) {
  193 + std::string msg = "Parameter \""+name+"\" not found in spatial "
  194 + "component of background model.";
  195 + throw GException::invalid_argument(G_ACCESS2, msg);
  196 + }
  197 +
  198 + // Return reference
  199 + return *(m_pars[index]);
  200 +}
  201 +
  202 +
  203 +/***********************************************************************//**
  204 + * @brief Returns reference to model parameter (const version)
  205 + *
  206 + * @param[in] name Parameter name.
  207 + *
  208 + * @exception GException::invalid_argument
  209 + * Parameter with specified name not found in container.
  210 + ***************************************************************************/
  211 +const GModelPar& GCTAModelSpatial::operator[](const std::string& name) const
  212 +{
  213 + // Get parameter index
  214 + int index = 0;
  215 + for (; index < size(); ++index) {
  216 + if (m_pars[index]->name() == name)
  217 + break;
  218 + }
  219 +
  220 + // Throw exception if parameter name was not found
  221 + if (index >= size()) {
  222 + std::string msg = "Parameter \""+name+"\" not found in spatial "
  223 + "component of background model.";
  224 + throw GException::invalid_argument(G_ACCESS2, msg);
  225 + }
  226 +
  227 + // Return reference
  228 + return *(m_pars[index]);
  229 +}
  230 +
  231 +
  232 +/*==========================================================================
  233 + = =
  234 + = Public methods =
  235 + = =
  236 + ==========================================================================*/
  237 +
  238 +/*==========================================================================
  239 + = =
  240 + = Private methods =
  241 + = =
  242 + ==========================================================================*/
  243 +
  244 +/***********************************************************************//**
  245 + * @brief Initialise class members
  246 + ***************************************************************************/
  247 +void GCTAModelSpatial::init_members(void)
  248 +{
  249 + // Initialise members
  250 + m_pars.clear();
  251 +
  252 + // Return
  253 + return;
  254 +}
  255 +
  256 +
  257 +/***********************************************************************//**
  258 + * @brief Copy class members
  259 + *
  260 + * @param[in] model Spatial acceptance model.
  261 + ***************************************************************************/
  262 +void GCTAModelSpatial::copy_members(const GCTAModelSpatial& model)
  263 +{
  264 + // Copy members
  265 + m_pars = model.m_pars;
  266 +
  267 + // Return
  268 + return;
  269 +}
  270 +
  271 +
  272 +/***********************************************************************//**
  273 + * @brief Delete class members
  274 + ***************************************************************************/
  275 +void GCTAModelSpatial::free_members(void)
  276 +{
  277 + // Return
  278 + return;
  279 +}
... ...
inst/cta/src/GCTAModelSpatialRegistry.cpp 0 → 100644
  View file @513f99b
  1 +/***************************************************************************
  2 + * GCTAModelSpatialRegistry.cpp - Spatial model registry class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2018 by Juergen Knoedlseder *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file GCTAModelSpatialRegistry.cpp
  23 + * @brief Spatial model registry class implemenation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GModelPar.hpp"
  32 +#include "GCTAModelSpatial.hpp"
  33 +#include "GCTAModelSpatialRegistry.hpp"
  34 +#include "GException.hpp"
  35 +#include "GTools.hpp"
  36 +
  37 +/* __ Method name definitions ____________________________________________ */
  38 +#define G_ALLOC "GCTAModelSpatialRegistry::alloc(GXmlElement&)"
  39 +#define G_NAME "GCTAModelSpatialRegistry::name(int&)"
  40 +
  41 +/* __ Macros _____________________________________________________________ */
  42 +
  43 +/* __ Coding definitions _________________________________________________ */
  44 +
  45 +/* __ Debug definitions __________________________________________________ */
  46 +//#define G_DEBUG_REGISTRY //!< Dump registry
  47 +
  48 +
  49 +/*==========================================================================
  50 + = =
  51 + = Constructors/destructors =
  52 + = =
  53 + ==========================================================================*/
  54 +
  55 +/***********************************************************************//**
  56 + * @brief Void constructor
  57 + ***************************************************************************/
  58 +GCTAModelSpatialRegistry::GCTAModelSpatialRegistry(void)
  59 +{
  60 + // Initialise private members for clean destruction
  61 + init_members();
  62 +
  63 + // Debug option: Show actual registry
  64 + #if defined(G_DEBUG_REGISTRY)
  65 + std::cout << "GCTAModelSpatialRegistry(void): ";
  66 + for (int i = 0; i < size(); ++i) {
  67 + std::cout << "\"" << models()[i]->type() << "\" ";
  68 + }
  69 + std::cout << std::endl;
  70 + #endif
  71 +
  72 + // Return
  73 + return;
  74 +}
  75 +
  76 +
  77 +/***********************************************************************//**
  78 + * @brief Model constructor
  79 + *
  80 + * @param[in] model Model.
  81 + *
  82 + * Construct registry by adding a model to the registry. This is the standard
  83 + * constructor that is used to register a new model to GammaLib.
  84 + ***************************************************************************/
  85 +GCTAModelSpatialRegistry::GCTAModelSpatialRegistry(const GCTAModelSpatial* model)
  86 +{
  87 + // Initialise private members for clean destruction
  88 + init_members();
  89 +
  90 + // Debug option: Notify new registry
  91 + #if defined(G_DEBUG_REGISTRY)
  92 + std::cout << "GCTAModelSpatialRegistry(const GCTAModelSpatial*): ";
  93 + std::cout << "add \"" << model->type() << "\" to registry." << std::endl;
  94 + #endif
  95 +
  96 + // Allocate new registry
  97 + const GCTAModelSpatial** new_models = new const GCTAModelSpatial*[size()+1];
  98 +
  99 + // Save old registry
  100 + for (int i = 0; i < size(); ++i) {
  101 + new_models[i] = models()[i];
  102 + }
  103 +
  104 + // Add new model to registry
  105 + new_models[size()] = model;
  106 +
  107 + // Set pointers on new registry
  108 + models().assign(new_models);
  109 +
  110 + // Increment number of models in registry
  111 + number()++;
  112 +
  113 + // Debug option: Show actual registry
  114 + #if defined(G_DEBUG_REGISTRY)
  115 + std::cout << "GCTAModelSpatialRegistry(const GCTAModelSpatial*): ";
  116 + for (int i = 0; i < size(); ++i) {
  117 + std::cout << "\"" << models()[i]->type() << "\" ";
  118 + }
  119 + std::cout << std::endl;
  120 + #endif
  121 +
  122 + // Return
  123 + return;
  124 +}
  125 +
  126 +
  127 +/***********************************************************************//**
  128 + * @brief Copy constructor
  129 + *
  130 + * @param[in] registry Registry.
  131 + ***************************************************************************/
  132 +GCTAModelSpatialRegistry::GCTAModelSpatialRegistry(const GCTAModelSpatialRegistry& registry)
  133 +{
  134 + // Initialise private members
  135 + init_members();
  136 +
  137 + // Copy members
  138 + copy_members(registry);
  139 +
  140 + // Return
  141 + return;
  142 +}
  143 +
  144 +
  145 +/***********************************************************************//**
  146 + * @brief Destructor
  147 + ***************************************************************************/
  148 +GCTAModelSpatialRegistry::~GCTAModelSpatialRegistry(void)
  149 +{
  150 + // Free members
  151 + free_members();
  152 +
  153 + // Return
  154 + return;
  155 +}
  156 +
  157 +
  158 +/*==========================================================================
  159 + = =
  160 + = Operators =
  161 + = =
  162 + ==========================================================================*/
  163 +
  164 +/***********************************************************************//**
  165 + * @brief Assignment operator
  166 + *
  167 + * @param[in] registry Registry.
  168 + * @return Reference to registry.
  169 + ***************************************************************************/
  170 +GCTAModelSpatialRegistry& GCTAModelSpatialRegistry::operator=(const GCTAModelSpatialRegistry& registry)
  171 +{
  172 + // Execute only if object is not identical
  173 + if (this != &registry) {
  174 +
  175 + // Free members
  176 + free_members();
  177 +
  178 + // Initialise private members
  179 + init_members();
  180 +
  181 + // Copy members
  182 + copy_members(registry);
  183 +
  184 + } // endif: object was not identical
  185 +
  186 + // Return
  187 + return *this;
  188 +}
  189 +
  190 +
  191 +/*==========================================================================
  192 + = =
  193 + = Public methods =
  194 + = =
  195 + ==========================================================================*/
  196 +
  197 +/***********************************************************************//**
  198 + * @brief Allocate spatial model that is found in XML element
  199 + *
  200 + * @param[in] xml XML element.
  201 + * @return Pointer to spatial model.
  202 + *
  203 + * @exception GException::invalid_value
  204 + * No appropriate spatial model found in XML element.
  205 + *
  206 + * Returns a pointer to a spatial model instance that corresponds to the
  207 + * type found in an XML element. If no appropriate model is found the method
  208 + * will throw an exception.
  209 + ***************************************************************************/
  210 +GCTAModelSpatial* GCTAModelSpatialRegistry::alloc(const GXmlElement& xml) const
  211 +{
  212 + // Initialise spatial model
  213 + GCTAModelSpatial* model = NULL;
  214 +
  215 + // Search for model type in registry
  216 + for (int i = 0; i < size(); ++i) {
  217 + if (models()[i]->type() == xml.attribute("type")) {
  218 + model = models()[i]->clone();
  219 + break;
  220 + }
  221 + }
  222 +
  223 + // If no model has been found then throw an exception
  224 + if (model == NULL) {
  225 + std::string msg = "Spatial model of type \""+xml.attribute("type")+
  226 + "\" not found in registry. Possible spatial model "
  227 + "types are:";
  228 + for (int i = 0; i < size(); ++i) {
  229 + msg += " \""+models()[i]->type()+"\"";
  230 + }
  231 + msg += ".";
  232 + throw GException::invalid_value(G_ALLOC, msg);
  233 + }
  234 +
  235 + // Read model from XML element
  236 + model->read(xml);
  237 +
  238 + // Return spectral model
  239 + return model;
  240 +}
  241 +
  242 +
  243 +/***********************************************************************//**
  244 + * @brief Returns model name
  245 + *
  246 + * @param[in] index Model index [0,...,size()-1].
  247 + * @return Model name.
  248 + *
  249 + * @exception GException::out_of_range
  250 + * Model index is out of range.
  251 + ***************************************************************************/
  252 +std::string GCTAModelSpatialRegistry::name(const int& index) const
  253 +{
  254 + // Compile option: raise exception if index is out of range
  255 + #if defined(G_RANGE_CHECK)
  256 + if (index < 0 || index >= size()) {
  257 + throw GException::out_of_range(G_NAME, "Spatial model index", index,
  258 + size());
  259 + }
  260 + #endif
  261 +
  262 + // Return name
  263 + return (models()[index]->type());
  264 +}
  265 +
  266 +
  267 +/***********************************************************************//**
  268 + * @brief Print registry information
  269 + *
  270 + * @param[in] chatter Chattiness (defaults to NORMAL).
  271 + * @return Registry content.
  272 + ***************************************************************************/
  273 +std::string GCTAModelSpatialRegistry::print(const GChatter& chatter) const
  274 +{
  275 + // Initialise result string
  276 + std::string result;
  277 +
  278 + // Continue only if chatter is not silent
  279 + if (chatter != SILENT) {
  280 +
  281 + // Append header
  282 + result.append("=== GCTAModelSpatialRegistry ===");
  283 +
  284 + // Append information
  285 + result.append("\n"+gammalib::parformat("Number of models"));
  286 + result.append(gammalib::str(size()));
  287 +
  288 + // NORMAL: Append models
  289 + if (chatter >= NORMAL) {
  290 + for (int i = 0; i < size(); ++i) {
  291 + result.append("\n"+gammalib::parformat(models()[i]->type()));
  292 + for (int k = 0; k < models()[i]->size(); ++k) {
  293 + if (k > 0) {
  294 + result.append(", ");
  295 + }
  296 + result.append("\"");
  297 + result.append((*(models()[i]))[k].name());
  298 + result.append("\"");
  299 + }
  300 + }
  301 + }
  302 +
  303 + } // endif: chatter was not silent
  304 +
  305 + // Return result
  306 + return result;
  307 +}
  308 +
  309 +/*==========================================================================
  310 + = =
  311 + = Private methods =
  312 + = =
  313 + ==========================================================================*/
  314 +
  315 +/***********************************************************************//**
  316 + * @brief Initialise class members
  317 + ***************************************************************************/
  318 +void GCTAModelSpatialRegistry::init_members(void)
  319 +{
  320 + // Return
  321 + return;
  322 +}
  323 +
  324 +
  325 +/***********************************************************************//**
  326 + * @brief Copy class members
  327 + *
  328 + * @param[in] registry Registry.
  329 + ***************************************************************************/
  330 +void GCTAModelSpatialRegistry::copy_members(const GCTAModelSpatialRegistry& registry)
  331 +{
  332 + // Return
  333 + return;
  334 +}
  335 +
  336 +
  337 +/***********************************************************************//**
  338 + * @brief Delete class members
  339 + ***************************************************************************/
  340 +void GCTAModelSpatialRegistry::free_members(void)
  341 +{
  342 + // Return
  343 + return;
  344 +}
... ...
inst/cta/test/data/cta_model_bgd_rad_gauss.xml 0 → 100644
  View file @513f99b
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="My model" type="CTABackground" instrument="CTA">
  4 + <spectrum type="PowerLaw">
  5 + <parameter name="Prefactor" scale="1e-6" value="61.8" min="0.0" max="1000.0" free="1"/>
  6 + <parameter name="Index" scale="-1" value="1.85" min="0.0" max="+5.0" free="1"/>
  7 + <parameter name="PivotEnergy" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
  8 + </spectrum>
  9 + <spatialModel type="Gaussian">
  10 + <parameter name="Sigma" scale="1.0" value="3.0" min="0.01" max="10.0" free="1"/>
  11 + </spatialModel>
  12 + </source>
  13 +</source_library>
... ...
inst/cta/test/test_CTA.cpp
  View file @513f99b
... ... @@ -47,29 +47,30 @@
47 47 /* __ Globals ____________________________________________________________ */
48 48  
49 49 /* __ Constants __________________________________________________________ */
50   -const std::string datadir = std::getenv("TEST_CTA_DATA");
51   -const std::string caldbdir = datadir + "/../caldb";
52   -const std::string cta_caldb = datadir + "/../../caldb";
53   -const std::string cta_irf = "cta_dummy_irf";
54   -const std::string cta_events = datadir+"/crab_events.fits";
55   -const std::string cta_events_gti = datadir+"/crab_events_gti.fits[EVENTS2]";
56   -const std::string cta_cntmap = datadir+"/crab_cntmap.fits";
57   -const std::string cta_bin_xml = datadir+"/obs_binned.xml";
58   -const std::string cta_unbin_xml = datadir+"/obs_unbinned.xml";
59   -const std::string cta_model_xml = datadir+"/crab.xml";
60   -const std::string cta_rsp_xml = datadir+"/rsp_models.xml";
61   -const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml";
62   -const std::string cta_irf_bgd_xml = datadir+"/cta_model_irf_bgd.xml";
63   -const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml";
64   -const std::string cta_caldb_king = cta_caldb+"/data/cta/e/bcf/IFAE20120510_50h_King";
65   -const std::string cta_irf_king = "irf_file.fits";
66   -const std::string cta_psf_table = caldbdir+"/psf_table.fits[PSF_2D_TABLE]";
67   -const std::string cta_perf_table = caldbdir+"/cta_dummy_irf.dat";
68   -const std::string cta_edisp_rmf = caldbdir+"/dc1/rmf.fits";
69   -const std::string cta_edisp_2D = caldbdir+"/edisp_matrix.fits";
70   -const std::string cta_bgd_3D = caldbdir+"/edisp_matrix.fits";
71   -const std::string cta_modbck_fit = datadir+"/bg_test.fits";
72   -const std::string cta_point_table = datadir+"/crab_pointing.fits";
  50 +const std::string datadir = std::getenv("TEST_CTA_DATA");
  51 +const std::string caldbdir = datadir + "/../caldb";
  52 +const std::string cta_caldb = datadir + "/../../caldb";
  53 +const std::string cta_irf = "cta_dummy_irf";
  54 +const std::string cta_events = datadir+"/crab_events.fits";
  55 +const std::string cta_events_gti = datadir+"/crab_events_gti.fits[EVENTS2]";
  56 +const std::string cta_cntmap = datadir+"/crab_cntmap.fits";
  57 +const std::string cta_bin_xml = datadir+"/obs_binned.xml";
  58 +const std::string cta_unbin_xml = datadir+"/obs_unbinned.xml";
  59 +const std::string cta_model_xml = datadir+"/crab.xml";
  60 +const std::string cta_rsp_xml = datadir+"/rsp_models.xml";
  61 +const std::string cta_bgd_rad_gauss_xml = datadir+"/cta_model_bgd_rad_gauss.xml";
  62 +const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml";
  63 +const std::string cta_irf_bgd_xml = datadir+"/cta_model_irf_bgd.xml";
  64 +const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml";
  65 +const std::string cta_caldb_king = cta_caldb+"/data/cta/e/bcf/IFAE20120510_50h_King";
  66 +const std::string cta_irf_king = "irf_file.fits";
  67 +const std::string cta_psf_table = caldbdir+"/psf_table.fits[PSF_2D_TABLE]";
  68 +const std::string cta_perf_table = caldbdir+"/cta_dummy_irf.dat";
  69 +const std::string cta_edisp_rmf = caldbdir+"/dc1/rmf.fits";
  70 +const std::string cta_edisp_2D = caldbdir+"/edisp_matrix.fits";
  71 +const std::string cta_bgd_3D = caldbdir+"/edisp_matrix.fits";
  72 +const std::string cta_modbck_fit = datadir+"/bg_test.fits";
  73 +const std::string cta_point_table = datadir+"/crab_pointing.fits";
73 74  
74 75 /* __ Test files for stacked analysis (based on Prod2::South_0.5h) _______ */
75 76 const std::string cta_stacked_xml = datadir+"/stacked_obs.xml";
... ... @@ -172,6 +173,8 @@ void TestGCTAModel::set(void)
172 173 name("Test CTA models");
173 174  
174 175 // Append tests to test suite
  176 + append(static_cast<pfunction>(&TestGCTAModel::test_model_bgd),
  177 + "Test CTA background model");
175 178 append(static_cast<pfunction>(&TestGCTAModel::test_model_cube_bgd),
176 179 "Test CTA cube background model");
177 180 append(static_cast<pfunction>(&TestGCTAModel::test_model_irf_bgd),
... ... @@ -1453,6 +1456,78 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
1453 1456  
1454 1457  
1455 1458 /***********************************************************************//**
  1459 + * @brief Test CTA background
  1460 + ***************************************************************************/
  1461 +void TestGCTAModel::test_model_bgd(void)
  1462 +{
  1463 + // Test void constuctor
  1464 + test_try("Test void constuctor");
  1465 + try {
  1466 + GCTAModelBackground model;
  1467 + test_try_success();
  1468 + }
  1469 + catch (std::exception &e) {
  1470 + test_try_failure(e);
  1471 + }
  1472 +
  1473 + // Test XML constuctor
  1474 + test_try("Test XML constuctor");
  1475 + try {
  1476 + GXml xml(cta_bgd_rad_gauss_xml);
  1477 + const GXmlElement& lib = *xml.element("source_library", 0);
  1478 + const GXmlElement& src = *lib.element("source", 0);
  1479 + GCTAModelBackground model(src);
  1480 + test_value(model["Prefactor"].value(), 61.8e-6);
  1481 + test_value(model["Index"].value(), -1.85);
  1482 + test_value(model["PivotEnergy"].value(), 1.0e6);
  1483 + test_assert(model.is_constant(), "Model is expected to be constant.");
  1484 + test_try_success();
  1485 + }
  1486 + catch (std::exception &e) {
  1487 + test_try_failure(e);
  1488 + }
  1489 +
  1490 + // Test spectral constuctor
  1491 + /*
  1492 + test_try("Test spectral constuctor");
  1493 + try {
  1494 + GModelSpectralPlaw plaw(1.0, 0.0, GEnergy(1.0, "TeV"));
  1495 + GCTAModelBackground model(plaw);
  1496 + test_value(model["Prefactor"].value(), 1.0);
  1497 + test_value(model["Index"].value(), 0.0);
  1498 + test_value(model["PivotEnergy"].value(), 1.0e6);
  1499 + test_assert(model.is_constant(), "Model is expected to be constant.");
  1500 + test_try_success();
  1501 + }
  1502 + catch (std::exception &e) {
  1503 + test_try_failure(e);
  1504 + }
  1505 + */
  1506 +
  1507 + // Test XML loading of cube background
  1508 + GModels models(cta_bgd_rad_gauss_xml);
  1509 + GModel* model = models["My model"];
  1510 + test_value((*model)["Prefactor"].value(), 61.8e-6);
  1511 + test_value((*model)["Index"].value(), -1.85);
  1512 + test_value((*model)["PivotEnergy"].value(), 1.0e6);
  1513 + test_assert(model->is_constant(), "Model is expected to be constant.");
  1514 +
  1515 + // Test XML saving and reloading of cube background
  1516 + models.save("test.xml");
  1517 + models.load("test.xml");
  1518 + model = models["My model"];
  1519 + test_value((*model)["Prefactor"].value(), 61.8e-6);
  1520 + test_value((*model)["Index"].value(), -1.85);
  1521 + test_value((*model)["PivotEnergy"].value(), 1.0e6);
  1522 + test_assert(model->is_constant(), "Model is expected to be constant.");
  1523 +
  1524 + // Return
  1525 + return;
  1526 +
  1527 +}
  1528 +
  1529 +
  1530 +/***********************************************************************//**
1456 1531 * @brief Test CTA cube background
1457 1532 ***************************************************************************/
1458 1533 void TestGCTAModel::test_model_cube_bgd(void)
... ...
inst/cta/test/test_CTA.hpp
  View file @513f99b
... ... @@ -118,6 +118,7 @@ public:
118 118 virtual void set(void);
119 119 virtual TestGCTAModel* clone(void) const { return new TestGCTAModel(*this); }
120 120 virtual std::string classname(void) const { return "TestGCTAModel"; }
  121 + void test_model_bgd(void);
121 122 void test_model_cube_bgd(void);
122 123 void test_model_irf_bgd(void);
123 124 void test_model_aeff_bgd(void);
... ...