Karim Sabri / gammalib

Commit c4b0b2aa004161b201aad4bab746ae8d9d06aee2

Authored by Jürgen Knödlseder
1 parent 7da0f9cd

Implement energy-dependent Gaussian CTA spatial background model (#2868)

Inline Side-by-side
Showing 13 changed files with 912 additions and 6 deletions
ChangeLog
  View file @c4b0b2a
1   -2019-04-09
  1 +2019-04-13
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
5 5  
  6 + Implement energy-dependent Gaussian CTA spatial background model (#2868)
6 7 Take instrument scale factors in On/Off analysis into account (#2867)
7 8 GWcs::wcs_set() method now sets radesys and equinox members (#2866)
8 9 Add GCTAPointing::is_valid() method and remove table methods
... ...
1 1 New Features and Important Changes in GammaLib 1.6.0
2 2  
3   -9 April 2019
  3 +13 April 2019
4 4  
5 5  
6 6 1. Introduction
... ... @@ -20,6 +20,7 @@ The following classes have been added:
20 20 - GCTAModelSpatial
21 21 - GCTAModelSpatialMultiplicative
22 22 - GCTAModelSpatialRegistry
  23 +- GCTAModelSpatialGaussSpectrum
23 24  
24 25 The following classes have been removed:
25 26 - GCTACubeSourceDiffuse
... ... @@ -331,6 +332,9 @@ unit tests for that response functions were added.
331 332  
332 333 20. CTA interface
333 334 -----------------
  335 +Add the GCTAModelSpatialGaussSpectrum class that implements an energy-dependent
  336 +Gaussian spatial background model (#2868).
  337 +
334 338 Instrument scale factors are now also taken in an On/Off analysis into account.
335 339 The scale factors were added to the computation in GCTAOnOffObservation::N_gamma()
336 340 (#2867).
... ...
doc/source/admin/release_history/1.6.rst
  View file @c4b0b2a
... ... @@ -55,6 +55,8 @@ Bug fixes
55 55 Improvements
56 56 ------------
57 57  
  58 +* [`2868 <https://cta-redmine.irap.omp.eu/issues/2868>`_] -
  59 + Implement energy-dependent Gaussian CTA spatial background model
58 60 * [`2867 <https://cta-redmine.irap.omp.eu/issues/2867>`_] -
59 61 Take instrument scale factors in On/Off analysis into account
60 62 * [`2850 <https://cta-redmine.irap.omp.eu/issues/2850>`_] -
... ...
inst/cta/Makefile.am
  View file @c4b0b2a
1 1 #############################################################################
2 2 # Makefile.am for cta module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2018 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2019 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -75,6 +75,7 @@ sources = src/GCTASupport.cpp \
75 75 src/GCTAPointing.cpp \
76 76 src/GCTAModelBackground.cpp \
77 77 src/GCTAModelSpatial.cpp \
  78 + src/GCTAModelSpatialGaussSpectrum.cpp \
78 79 src/GCTAModelSpatialGradient.cpp \
79 80 src/GCTAModelSpatialMultiplicative.cpp \
80 81 src/GCTAModelSpatialRegistry.cpp \
... ... @@ -129,6 +130,7 @@ pkginclude_HEADERS = include/GCTATypemaps.hpp \
129 130 include/GCTACubeSource.hpp \
130 131 include/GCTAModelBackground.hpp \
131 132 include/GCTAModelSpatial.hpp \
  133 + include/GCTAModelSpatialGaussSpectrum.hpp \
132 134 include/GCTAModelSpatialGradient.hpp \
133 135 include/GCTAModelSpatialMultiplicative.hpp \
134 136 include/GCTAModelSpatialRegistry.hpp \
... ...
inst/cta/include/GCTALib.hpp
  View file @c4b0b2a
1 1 /***************************************************************************
2 2 * GCTALib.hpp - CTA Support Header files *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2019 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -73,6 +73,7 @@
73 73 #include "GCTACubeSource.hpp"
74 74 #include "GCTAModelBackground.hpp"
75 75 #include "GCTAModelSpatial.hpp"
  76 +#include "GCTAModelSpatialGaussSpectrum.hpp"
76 77 #include "GCTAModelSpatialGradient.hpp"
77 78 #include "GCTAModelSpatialMultiplicative.hpp"
78 79 #include "GCTAModelSpatialRegistry.hpp"
... ...
inst/cta/include/GCTAModelSpatialGaussSpectrum.hpp 0 → 100644
  View file @c4b0b2a
  1 +/***************************************************************************
  2 + * GCTAModelSpatialGaussSpectrum.hpp - Spatial energy dependent Gaussian *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2019 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 GCTAModelSpatialGaussSpectrum.hpp
  23 + * @brief Spatial energy dependent Gaussian interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCTAMODELSPATIALGAUSSSPECTRUM_HPP
  28 +#define GCTAMODELSPATIALGAUSSSPECTRUM_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GModelSpectral.hpp"
  33 +#include "GCTAModelSpatial.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GXmlElement;
  37 +class GCTAInstDir;
  38 +class GEnergy;
  39 +class GTime;
  40 +class GCTAObservation;
  41 +
  42 +
  43 +/***********************************************************************//**
  44 + * @class GCTAModelSpatialGaussSpectrum
  45 + *
  46 + * @brief Spatial energy dependent Gaussian model class
  47 + ***************************************************************************/
  48 +class GCTAModelSpatialGaussSpectrum : public GCTAModelSpatial {
  49 +
  50 +public:
  51 + // Constructors and destructors
  52 + GCTAModelSpatialGaussSpectrum(void);
  53 + explicit GCTAModelSpatialGaussSpectrum(const double& sigma);
  54 + explicit GCTAModelSpatialGaussSpectrum(const GModelSpectral& sigma);
  55 + explicit GCTAModelSpatialGaussSpectrum(const GXmlElement& xml);
  56 + GCTAModelSpatialGaussSpectrum(const GCTAModelSpatialGaussSpectrum& model);
  57 + virtual ~GCTAModelSpatialGaussSpectrum(void);
  58 +
  59 + // Operators
  60 + virtual GCTAModelSpatialGaussSpectrum& operator=(const GCTAModelSpatialGaussSpectrum& model);
  61 +
  62 + // Implemented pure virtual methods
  63 + virtual void clear(void);
  64 + virtual GCTAModelSpatialGaussSpectrum* clone(void) const;
  65 + virtual std::string classname(void) const;
  66 + virtual std::string type(void) const;
  67 + virtual double eval(const GCTAInstDir& dir,
  68 + const GEnergy& energy,
  69 + const GTime& time,
  70 + const bool& gradients = false) const;
  71 + virtual double mc_max_value(const GCTAObservation& obs) const;
  72 + virtual void read(const GXmlElement& xml);
  73 + virtual void write(GXmlElement& xml) const;
  74 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  75 +
  76 + // Other methods
  77 + const GModelSpectral* sigma(void) const;
  78 + void sigma(const double& sigma);
  79 + void sigma(const GModelSpectral& sigma);
  80 +
  81 +protected:
  82 + // Protected methods
  83 + void init_members(void);
  84 + void copy_members(const GCTAModelSpatialGaussSpectrum& model);
  85 + void free_members(void);
  86 + void set_pointers(void);
  87 +
  88 +
  89 + // Protected members
  90 + GModelSpectral* m_sigma; // Pointer to sigma spectrum
  91 +};
  92 +
  93 +
  94 +/***********************************************************************//**
  95 + * @brief Return class name
  96 + *
  97 + * @return String containing the class name ("GCTAModelSpatialGaussSpectrum").
  98 + ***************************************************************************/
  99 +inline
  100 +std::string GCTAModelSpatialGaussSpectrum::classname(void) const
  101 +{
  102 + return ("GCTAModelSpatialGaussSpectrum");
  103 +}
  104 +
  105 +
  106 +/***********************************************************************//**
  107 + * @brief Return model type
  108 + *
  109 + * @return Model type "EnergyDependentGaussian".
  110 + ***************************************************************************/
  111 +inline
  112 +std::string GCTAModelSpatialGaussSpectrum::type(void) const
  113 +{
  114 + return ("EnergyDependentGaussian");
  115 +}
  116 +
  117 +
  118 +/***********************************************************************//**
  119 + * @brief Return maximum function value for Monte Carlo simulations
  120 + *
  121 + * @param[in] obs CTA Observation.
  122 + * @return Maximum function value for Monte Carlo simulations.
  123 + *
  124 + * This method always returns 1.
  125 + ***************************************************************************/
  126 +inline
  127 +double GCTAModelSpatialGaussSpectrum::mc_max_value(const GCTAObservation& obs) const
  128 +{
  129 + return 1.0;
  130 +}
  131 +
  132 +
  133 +/***********************************************************************//**
  134 + * @brief Return pointer to sigma spectrum
  135 + *
  136 + * @return Pointer to sigma spectrum.
  137 + ***************************************************************************/
  138 +inline
  139 +const GModelSpectral* GCTAModelSpatialGaussSpectrum::sigma(void) const
  140 +{
  141 + return (m_sigma);
  142 +}
  143 +
  144 +#endif /* GCTAMODELSPATIALGAUSSSPECTRUM_HPP */
... ...
inst/cta/pyext/GCTAModelSpatialGaussSpectrum.i 0 → 100644
  View file @c4b0b2a
  1 +/***************************************************************************
  2 + * GCTAModelSpatialGaussSpectrum.i - Spatial energy dependent Gaussian *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2019 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 GCTAModelSpatialGaussSpectrum.i
  23 + * @brief Spatial energy dependent Gaussian interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCTAModelSpatialGaussSpectrum.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCTAModelSpatialGaussSpectrum
  34 + *
  35 + * @brief Spatial energy dependent Gaussian model class
  36 + ***************************************************************************/
  37 +class GCTAModelSpatialGaussSpectrum : public GCTAModelSpatial {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + GCTAModelSpatialGaussSpectrum(void);
  42 + explicit GCTAModelSpatialGaussSpectrum(const double& sigma);
  43 + explicit GCTAModelSpatialGaussSpectrum(const GModelSpectral& sigma);
  44 + explicit GCTAModelSpatialGaussSpectrum(const GXmlElement& xml);
  45 + GCTAModelSpatialGaussSpectrum(const GCTAModelSpatialGaussSpectrum& model);
  46 + virtual ~GCTAModelSpatialGaussSpectrum(void);
  47 +
  48 + // Operators
  49 + virtual GCTAModelSpatialGaussSpectrum& operator=(const GCTAModelSpatialGaussSpectrum& model);
  50 +
  51 + // Implemented pure virtual methods
  52 + virtual void clear(void);
  53 + virtual GCTAModelSpatialGaussSpectrum* clone(void) const;
  54 + virtual std::string classname(void) const;
  55 + virtual std::string type(void) const;
  56 + virtual double eval(const GCTAInstDir& dir,
  57 + const GEnergy& energy,
  58 + const GTime& time,
  59 + const bool& gradients = false) const;
  60 + virtual double mc_max_value(const GCTAObservation& obs) const;
  61 + virtual void read(const GXmlElement& xml);
  62 + virtual void write(GXmlElement& xml) const;
  63 +
  64 + // Other methods
  65 + const GModelSpectral* sigma(void) const;
  66 + void sigma(const double& sigma);
  67 + void sigma(const GModelSpectral& sigma);
  68 +};
  69 +
  70 +
  71 +/***********************************************************************//**
  72 + * @brief GCTAModelSpatialGaussSpectrum class extension
  73 + ***************************************************************************/
  74 +%extend GCTAModelSpatialGaussSpectrum {
  75 + GCTAModelSpatialGaussSpectrum copy() {
  76 + return (*self);
  77 + }
  78 +%pythoncode {
  79 + def __getstate__(self):
  80 + xml = gammalib.GXmlElement()
  81 + self.write(xml)
  82 + state = (xml,)
  83 + return state
  84 + def __setstate__(self, state):
  85 + if state[0].elements('parameter') == 0:
  86 + self.__init__()
  87 + else:
  88 + self.__init__(state[0])
  89 +}
  90 +};
... ...
inst/cta/pyext/cta.i
  View file @c4b0b2a
1 1 /***************************************************************************
2 2 * cta.i - Cherenkov Telescope Array module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2018 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2019 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -308,6 +308,7 @@
308 308 %include "GCTAModelIrfBackground.i"
309 309 %include "GCTAModelAeffBackground.i"
310 310 %include "GCTAModelSpatial.i"
  311 +%include "GCTAModelSpatialGaussSpectrum.i"
311 312 %include "GCTAModelSpatialGradient.i"
312 313 %include "GCTAModelSpatialMultiplicative.i"
313 314 %include "GCTAModelSpatialRegistry.i"
... ...
inst/cta/src/GCTAModelRadialGauss.cpp
  View file @c4b0b2a
... ... @@ -533,7 +533,7 @@ void GCTAModelRadialGauss::write(GXmlElement&amp; xml) const
533 533 /***********************************************************************//**
534 534 * @brief Print point source information
535 535 *
536   - * @param[in] chatter Chattiness (defaults to NORMAL).
  536 + * @param[in] chatter Chattiness.
537 537 * @return String containing point source information.
538 538 ***************************************************************************/
539 539 std::string GCTAModelRadialGauss::print(const GChatter& chatter) const
... ...
inst/cta/src/GCTAModelSpatialGaussSpectrum.cpp 0 → 100644
  View file @c4b0b2a
  1 +/***************************************************************************
  2 + * GCTAModelSpatialGaussSpectrum.cpp - Spatial energy dependent Gaussian *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2019 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 GCTAModelSpatialGaussSpectrum.cpp
  23 + * @brief Spatial energy dependent Gaussian interface implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cmath>
  32 +#include "GTools.hpp"
  33 +#include "GMath.hpp"
  34 +#include "GCTAInstDir.hpp"
  35 +#include "GEnergy.hpp"
  36 +#include "GTime.hpp"
  37 +#include "GCTAModelSpatialGaussSpectrum.hpp"
  38 +#include "GCTAModelSpatialRegistry.hpp"
  39 +#include "GModelSpectralRegistry.hpp"
  40 +#include "GModelSpectralConst.hpp"
  41 +
  42 +/* __ Constants __________________________________________________________ */
  43 +
  44 +/* __ Globals ____________________________________________________________ */
  45 +const GCTAModelSpatialGaussSpectrum g_cta_spatial_gauss_spec_seed;
  46 +const GCTAModelSpatialRegistry g_cta_spatial_gauss_spec_registry(&g_cta_spatial_gauss_spec_seed);
  47 +
  48 +/* __ Method name definitions ____________________________________________ */
  49 +#define G_READ "GCTAModelSpatialGaussSpectrum::read(GXmlElement&)"
  50 +#define G_WRITE "GCTAModelSpatialGaussSpectrum::write(GXmlElement&)"
  51 +
  52 +/* __ Macros _____________________________________________________________ */
  53 +
  54 +/* __ Coding definitions _________________________________________________ */
  55 +
  56 +/* __ Debug definitions __________________________________________________ */
  57 +
  58 +
  59 +/*==========================================================================
  60 + = =
  61 + = Constructors/destructors =
  62 + = =
  63 + ==========================================================================*/
  64 +
  65 +/***********************************************************************//**
  66 + * @brief Void constructor
  67 + ***************************************************************************/
  68 +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(void) :
  69 + GCTAModelSpatial()
  70 +{
  71 + // Initialise members
  72 + init_members();
  73 +
  74 + // Return
  75 + return;
  76 +}
  77 +
  78 +
  79 +/***********************************************************************//**
  80 + * @brief Energy-independent Gaussian constructor
  81 + *
  82 + * @param[in] sigma Gaussian width (degrees\f$^2\f$).
  83 + ***************************************************************************/
  84 +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(const double& sigma) :
  85 + GCTAModelSpatial()
  86 +{
  87 + // Initialise members
  88 + init_members();
  89 +
  90 + // Assign sigma
  91 + this->sigma(sigma);
  92 +
  93 + // Return
  94 + return;
  95 +}
  96 +
  97 +
  98 +/***********************************************************************//**
  99 + * @brief Gaussian spectrum constructor
  100 + *
  101 + * @param[in] sigma Spectrum defining the energy-deependent Gaussian width
  102 + * (degrees\f$^2\f$).
  103 + ***************************************************************************/
  104 +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(const GModelSpectral& sigma) :
  105 + GCTAModelSpatial()
  106 +{
  107 + // Initialise members
  108 + init_members();
  109 +
  110 + // Assign sigma spectrum
  111 + this->sigma(sigma);
  112 +
  113 + // Return
  114 + return;
  115 +}
  116 +
  117 +
  118 +/***********************************************************************//**
  119 + * @brief XML constructor
  120 + *
  121 + * @param[in] xml XML element.
  122 + *
  123 + * Creates instance of a spatial energy-dependent Gaussian model by
  124 + * extracting information from an XML element. See
  125 + * GCTAModelSpatialGaussSpectrum::read() for more information about the
  126 + * expected structure of the XML element.
  127 + ***************************************************************************/
  128 +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(const GXmlElement& xml) :
  129 + GCTAModelSpatial()
  130 +{
  131 + // Initialise members
  132 + init_members();
  133 +
  134 + // Read information from XML element
  135 + read(xml);
  136 +
  137 + // Return
  138 + return;
  139 +}
  140 +
  141 +
  142 +/***********************************************************************//**
  143 + * @brief Copy constructor
  144 + *
  145 + * @param[in] model Energy-dependent Gaussian model.
  146 + ***************************************************************************/
  147 +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(const GCTAModelSpatialGaussSpectrum& model) :
  148 + GCTAModelSpatial(model)
  149 +{
  150 + // Initialise members
  151 + init_members();
  152 +
  153 + // Copy members
  154 + copy_members(model);
  155 +
  156 + // Return
  157 + return;
  158 +}
  159 +
  160 +
  161 +/***********************************************************************//**
  162 + * @brief Destructor
  163 + ***************************************************************************/
  164 +GCTAModelSpatialGaussSpectrum::~GCTAModelSpatialGaussSpectrum(void)
  165 +{
  166 + // Free members
  167 + free_members();
  168 +
  169 + // Return
  170 + return;
  171 +}
  172 +
  173 +
  174 +/*==========================================================================
  175 + = =
  176 + = Operators =
  177 + = =
  178 + ==========================================================================*/
  179 +
  180 +/***********************************************************************//**
  181 + * @brief Assignment operator
  182 + *
  183 + * @param[in] model Energy-dependent Gaussian model.
  184 + ***************************************************************************/
  185 +GCTAModelSpatialGaussSpectrum& GCTAModelSpatialGaussSpectrum::operator=(const GCTAModelSpatialGaussSpectrum& model)
  186 +{
  187 + // Execute only if object is not identical
  188 + if (this != &model) {
  189 +
  190 + // Copy base class members
  191 + this->GCTAModelSpatial::operator=(model);
  192 +
  193 + // Free members
  194 + free_members();
  195 +
  196 + // Initialise members
  197 + init_members();
  198 +
  199 + // Copy members
  200 + copy_members(model);
  201 +
  202 + } // endif: object was not identical
  203 +
  204 + // Return
  205 + return *this;
  206 +}
  207 +
  208 +
  209 +/*==========================================================================
  210 + = =
  211 + = Public methods =
  212 + = =
  213 + ==========================================================================*/
  214 +
  215 +/***********************************************************************//**
  216 + * @brief Clear instance
  217 + ***************************************************************************/
  218 +void GCTAModelSpatialGaussSpectrum::clear(void)
  219 +{
  220 + // Free class members (base and derived classes, derived class first)
  221 + free_members();
  222 + this->GCTAModelSpatial::free_members();
  223 +
  224 + // Initialise members
  225 + this->GCTAModelSpatial::init_members();
  226 + init_members();
  227 +
  228 + // Return
  229 + return;
  230 +}
  231 +
  232 +
  233 +/***********************************************************************//**
  234 + * @brief Clone instance
  235 + ***************************************************************************/
  236 +GCTAModelSpatialGaussSpectrum* GCTAModelSpatialGaussSpectrum::clone(void) const
  237 +{
  238 + return new GCTAModelSpatialGaussSpectrum(*this);
  239 +}
  240 +
  241 +
  242 +/***********************************************************************//**
  243 + * @brief Evaluate function
  244 + *
  245 + * @param[in] dir Event direction.
  246 + * @param[in] energy Event energy.
  247 + * @param[in] time Event time.
  248 + * @param[in] gradients Compute gradients?
  249 + * @return Function value
  250 + *
  251 + * Evaluates the energy-dependent Gaussian model for a given event direction.
  252 + * The Gaussian model is defined by
  253 + *
  254 + * \f[f(\theta,E) = \exp \left(-\frac{1}{2}
  255 + * \left( \frac{\theta^2}{\sigma(E)} \right)^2 \right)\f]
  256 + *
  257 + * where
  258 + * \f$\theta\f$ is the offset angle (in degrees), and
  259 + * \f$\sigma(E)\f$ is the energy-dependent Gaussian width (in degrees\f$^2\f$).
  260 + *
  261 + * If the @p gradients flag is true the method will also compute the partial
  262 + * derivatives of the parameters. The partial derivative of the Gaussian width
  263 + * is given by
  264 + *
  265 + * \f[\frac{df}{d\sigma_v} = f(\theta) \frac{\theta^4}{\sigma^3} \sigma_s\f]
  266 + *
  267 + * where
  268 + * \f$\sigma_v\f$ is the value part,
  269 + * \f$\sigma_s\f$ is the scaling part, and
  270 + * \f$\sigma = \sigma_v \sigma_s\f$.
  271 + *
  272 + * Note that this method implements a function which is unity for
  273 + * \f$\theta=0\f$.
  274 + ***************************************************************************/
  275 +double GCTAModelSpatialGaussSpectrum::eval(const GCTAInstDir& dir,
  276 + const GEnergy& energy,
  277 + const GTime& time,
  278 + const bool& gradients) const
  279 +{
  280 + // Initialise value
  281 + double value = 0.0;
  282 +
  283 + // Continue only if a sigma spectrum exists
  284 + if (m_sigma != NULL) {
  285 +
  286 + // Compute offset angle in degrees
  287 + double offset = dir.theta() * gammalib::rad2deg;
  288 +
  289 + // Compute energy dependent sigma value in degrees^2
  290 + double sigma = m_sigma->eval(energy, time, gradients);
  291 +
  292 + // Compute Gaussian value
  293 + double arg = offset * offset / sigma;
  294 + double arg2 = arg * arg;
  295 + value = std::exp(-0.5 * arg2);
  296 +
  297 + // Optionally compute partial derivatives
  298 + if (gradients) {
  299 +
  300 + // Compute partial derivatives of the sigma parameter
  301 + double g_sigma = value * arg2 / sigma;
  302 +
  303 + // Loop over model parameters
  304 + for (int i = 0; i < m_sigma->size(); ++i) {
  305 +
  306 + // Get reference to model parameter
  307 + GModelPar& par = m_sigma->operator[](i);
  308 +
  309 + // Scale parameter gradient
  310 + par.gradient(par.gradient()*g_sigma);
  311 +
  312 + } // endfor: loop over model parameters
  313 +
  314 + } // endif: computed partial derivatives
  315 +
  316 + // Compile option: Check for NaN/Inf
  317 + #if defined(G_NAN_CHECK)
  318 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  319 + std::cout << "*** ERROR: GCTAModelSpatialGaussSpectrum::eval";
  320 + std::cout << "(offset=" << offset << "): NaN/Inf encountered";
  321 + std::cout << " (value=" << value;
  322 + std::cout << ", energy=" << energy;
  323 + std::cout << ", sigma=" << sigma;
  324 + std::cout << ")" << std::endl;
  325 + }
  326 + #endif
  327 +
  328 + } // endif: sigma spectrum existed
  329 +
  330 + // Return value
  331 + return value;
  332 +}
  333 +
  334 +
  335 +/***********************************************************************//**
  336 + * @brief Read model from XML element
  337 + *
  338 + * @param[in] xml XML element.
  339 + *
  340 + * Read the energy-dependent Gaussian spatial model information from an XML
  341 + * element. The XML element needs to be of the following format:
  342 + *
  343 + * <spatialModel type="EnergyDependentGaussian">
  344 + * <sigma type="...">
  345 + * ...
  346 + * </sigma>
  347 + * </spatialModel>
  348 + *
  349 + * where any spectral model type can be provided for the @a sigma tag.
  350 + ***************************************************************************/
  351 +void GCTAModelSpatialGaussSpectrum::read(const GXmlElement& xml)
  352 +{
  353 + // Clear model
  354 + clear();
  355 +
  356 + // Get pointers on sigma model components
  357 + const GXmlElement* sigma = xml.element("sigma", 0);
  358 +
  359 + // Get sigma model
  360 + GModelSpectralRegistry registry;
  361 + m_sigma = registry.alloc(*sigma);
  362 +
  363 + // Set pointers
  364 + set_pointers();
  365 +
  366 + // Return
  367 + return;
  368 +}
  369 +
  370 +
  371 +/***********************************************************************//**
  372 + * @brief Write model into XML element
  373 + *
  374 + * @param[in] xml XML element.
  375 + *
  376 + * @exception GException::invalid_value
  377 + * Spatial model is not of valid type.
  378 + *
  379 + * Write the energy-dependent Gaussian spatial model information into an XML
  380 + * element. The XML element will be of the following format:
  381 + *
  382 + * <spatialModel type="EnergyDependentGaussian">
  383 + * <sigma type="...">
  384 + * ...
  385 + * </sigma>
  386 + * </spatialModel>
  387 + *
  388 + * where the spectral sigma model will be written into the @a sigma tag.
  389 + ***************************************************************************/
  390 +void GCTAModelSpatialGaussSpectrum::write(GXmlElement& xml) const
  391 +{
  392 + // Set model type
  393 + if (xml.attribute("type") == "") {
  394 + xml.attribute("type", type());
  395 + }
  396 +
  397 + // Verify model type
  398 + if (xml.attribute("type") != type()) {
  399 + std::string msg = "Spatial model \""+xml.attribute("type")+
  400 + "\" is not of type \""+type()+"\".";
  401 + throw GException::invalid_value(G_WRITE, msg);
  402 + }
  403 +
  404 + // Write sigma spectrum if it exists
  405 + if (m_sigma != NULL) {
  406 +
  407 + // Create new sigma spectrum node
  408 + xml.append(GXmlElement("sigma"));
  409 +
  410 + // Get new sigma spectrum node
  411 + GXmlElement* sigma = xml.element("sigma", xml.elements("sigma")-1);
  412 +
  413 + // Write sigma spectrum
  414 + m_sigma->write(*sigma);
  415 +
  416 + } // endif: sigma spectrum was not NULL
  417 +
  418 + // Return
  419 + return;
  420 +}
  421 +
  422 +
  423 +/***********************************************************************//**
  424 + * @brief Set sigma spectrum from value
  425 + *
  426 + * @param[in] sigma Constant sigma value.
  427 + ***************************************************************************/
  428 +void GCTAModelSpatialGaussSpectrum::sigma(const double& sigma)
  429 +{
  430 + // Clear sigma spectrum
  431 + clear();
  432 +
  433 + // Set constant spectrum
  434 + m_sigma = new GModelSpectralConst(sigma);
  435 +
  436 + // Set parameter pointers
  437 + set_pointers();
  438 +
  439 + // Return
  440 + return;
  441 +}
  442 +
  443 +
  444 +/***********************************************************************//**
  445 + * @brief Set sigma spectrum
  446 + *
  447 + * @param[in] sigma Sigma spectrum.
  448 + ***************************************************************************/
  449 +void GCTAModelSpatialGaussSpectrum::sigma(const GModelSpectral& sigma)
  450 +{
  451 + // Clear sigma spectrum
  452 + clear();
  453 +
  454 + // Set spectrum
  455 + m_sigma = sigma.clone();
  456 +
  457 + // Set parameter pointers
  458 + set_pointers();
  459 +
  460 + // Return
  461 + return;
  462 +}
  463 +
  464 +
  465 +/***********************************************************************//**
  466 + * @brief Print point source information
  467 + *
  468 + * @param[in] chatter Chattiness.
  469 + * @return String containing point source information.
  470 + ***************************************************************************/
  471 +std::string GCTAModelSpatialGaussSpectrum::print(const GChatter& chatter) const
  472 +{
  473 + // Initialise result string
  474 + std::string result;
  475 +
  476 + // Continue only if chatter is not silent
  477 + if (chatter != SILENT) {
  478 +
  479 + // Append header
  480 + result.append("=== GCTAModelSpatialGaussSpectrum ===");
  481 +
  482 + // Determine number of sigma parameters
  483 + int n_sigma = (m_sigma != NULL) ? m_sigma->size() : 0;
  484 +
  485 + // Append sigma spectrum type
  486 + result.append("\n"+gammalib::parformat("Sigma spectrum"));
  487 + if (n_sigma > 0) {
  488 + result.append("\""+m_sigma->type()+"\"");
  489 + }
  490 +
  491 + // Append sigma spectrum parameters
  492 + for (int i = 0; i < n_sigma; ++i) {
  493 + result.append("\n"+(m_sigma)[i].print());
  494 + }
  495 +
  496 + } // endif: chatter was not silent
  497 +
  498 + // Return result
  499 + return result;
  500 +}
  501 +
  502 +
  503 +/*==========================================================================
  504 + = =
  505 + = Private methods =
  506 + = =
  507 + ==========================================================================*/
  508 +
  509 +/***********************************************************************//**
  510 + * @brief Initialise class members
  511 + ***************************************************************************/
  512 +void GCTAModelSpatialGaussSpectrum::init_members(void)
  513 +{
  514 + // Initialise sigma spectrum pointer
  515 + m_sigma = NULL;
  516 +
  517 + // Clear parameter pointer(s)
  518 + m_pars.clear();
  519 +
  520 + // Return
  521 + return;
  522 +}
  523 +
  524 +
  525 +/***********************************************************************//**
  526 + * @brief Copy class members
  527 + *
  528 + * @param[in] model Energy-dependent Gaussian model.
  529 + ***************************************************************************/
  530 +void GCTAModelSpatialGaussSpectrum::copy_members(const GCTAModelSpatialGaussSpectrum& model)
  531 +{
  532 + // Clone sigma spectrum
  533 + if (model.m_sigma != NULL) {
  534 + m_sigma = model.m_sigma->clone();
  535 + }
  536 +
  537 + // Set parameter pointers
  538 + set_pointers();
  539 +
  540 + // Return
  541 + return;
  542 +}
  543 +
  544 +
  545 +/***********************************************************************//**
  546 + * @brief Delete class members
  547 + ***************************************************************************/
  548 +void GCTAModelSpatialGaussSpectrum::free_members(void)
  549 +{
  550 + // Delete sigma spectrum
  551 + if (m_sigma != NULL) delete m_sigma;
  552 +
  553 + // Signal free pointer
  554 + m_sigma = NULL;
  555 +
  556 + // Return
  557 + return;
  558 +}
  559 +
  560 +
  561 +/***********************************************************************//**
  562 + * @brief Set pointers
  563 + *
  564 + * Set pointers to all model parameters. The pointers are stored in a vector
  565 + * that is member of the GModelData base class.
  566 + ***************************************************************************/
  567 +void GCTAModelSpatialGaussSpectrum::set_pointers(void)
  568 +{
  569 + // Clear parameters
  570 + m_pars.clear();
  571 +
  572 + // Determine number of sigma parameters
  573 + int n_sigma = (m_sigma != NULL) ? m_sigma->size() : 0;
  574 +
  575 + // Gather sigma parameters if there are some
  576 + if (n_sigma > 0) {
  577 + for (int i = 0; i < n_sigma; ++i) {
  578 + m_pars.push_back(&((*m_sigma)[i]));
  579 + }
  580 + }
  581 +
  582 + // Return
  583 + return;
  584 +}
... ...
inst/cta/test/data/cta_model_bgd_gauss_e.xml 0 → 100644
  View file @c4b0b2a
  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="EnergyDependentGaussian">
  10 + <sigma type="NodeFunction">
  11 + <node>
  12 + <parameter name="Energy" value="1" scale="300000" free="0" />
  13 + <parameter name="Intensity" value="3.0" scale="1.0" min="1.5" max="100.0" free="1" />
  14 + </node>
  15 + <node>
  16 + <parameter name="Energy" value="1" scale="10000000" free="0" />
  17 + <parameter name="Intensity" value="4.0" scale="1.0" min="1.5" max="100.0" free="1" />
  18 + </node>
  19 + </sigma>
  20 + </spatialModel>
  21 + </source>
  22 +</source_library>
... ...
inst/cta/test/test_CTA.cpp
  View file @c4b0b2a
... ... @@ -59,6 +59,7 @@ const std::string cta_unbin_xml = datadir+&quot;/obs_unbinned.xml&quot;;
59 59 const std::string cta_model_xml = datadir+"/crab.xml";
60 60 const std::string cta_rsp_xml = datadir+"/rsp_models.xml";
61 61 const std::string cta_bgd_rad_gauss_xml = datadir+"/cta_model_bgd_rad_gauss.xml";
  62 +const std::string cta_bgd_gauss_e_xml = datadir+"/cta_model_bgd_gauss_e.xml";
62 63 const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml";
63 64 const std::string cta_irf_bgd_xml = datadir+"/cta_model_irf_bgd.xml";
64 65 const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml";
... ... @@ -176,6 +177,8 @@ void TestGCTAModel::set(void)
176 177 "Test CTA IRF background model");
177 178 append(static_cast<pfunction>(&TestGCTAModel::test_model_aeff_bgd),
178 179 "Test CTA Aeff background model");
  180 + append(static_cast<pfunction>(&TestGCTAModel::test_spatial_gauss_spectrum),
  181 + "Test GCTAModelSpatialGaussSpectrum class");
179 182  
180 183 // Return
181 184 return;
... ... @@ -2029,6 +2032,57 @@ void TestGCTAModel::test_model_aeff_bgd(void)
2029 2032  
2030 2033  
2031 2034 /***********************************************************************//**
  2035 + * @brief Test GCTAModelSpatialGaussSpectrum class
  2036 + ***************************************************************************/
  2037 +void TestGCTAModel::test_spatial_gauss_spectrum(void)
  2038 +{
  2039 + // Set instrument direction, energy and time
  2040 + GCTAInstDir dir(0.0, 0.0);
  2041 + GEnergy energy(1.0, "TeV");
  2042 + GTime time;
  2043 +
  2044 + // Test void constructor
  2045 + GCTAModelSpatialGaussSpectrum model1;
  2046 + test_assert(model1.sigma() == NULL, "Test sigma() for void constructor");
  2047 + test_value(model1.eval(dir, energy, time), 0.0, "Test eval() for void constructor");
  2048 +
  2049 + // Test value constructor
  2050 + GCTAModelSpatialGaussSpectrum model2(3.0);
  2051 + test_assert(model2.sigma() != NULL, "Test sigma() for value constructor");
  2052 + test_value(model2.eval(dir, energy, time), 1.0, "Test eval() for value constructor");
  2053 +
  2054 + // Test spectrum constructor
  2055 + GModelSpectralConst spectrum(3.0);
  2056 + GCTAModelSpatialGaussSpectrum model3(spectrum);
  2057 + test_assert(model3.sigma() != NULL, "Test sigma() for spectrum constructor");
  2058 + test_value(model3.eval(dir, energy, time), 1.0, "Test eval() for spectrum constructor");
  2059 +
  2060 + // Test XML constructor
  2061 + GXml xml(cta_bgd_gauss_e_xml);
  2062 + GXmlElement* element = xml.element("source_library > source[0] > spatialModel");
  2063 + GCTAModelSpatialGaussSpectrum model4(*element);
  2064 + test_assert(model4.sigma() != NULL, "Test sigma() for XML constructor");
  2065 + test_value(model4.eval(dir, energy, time), 1.0, "Test eval() for XML constructor");
  2066 +
  2067 + // Test copy constructor
  2068 + GCTAModelSpatialGaussSpectrum model5(model4);
  2069 + test_assert(model5.sigma() != NULL, "Test sigma() for copy constructor");
  2070 + test_value(model5.eval(dir, energy, time), 1.0, "Test eval() for copy constructor");
  2071 +
  2072 + // Test sigma() setter
  2073 + GCTAModelSpatialGaussSpectrum model6;
  2074 + test_assert(model6.sigma() == NULL, "Test sigma() setter");
  2075 + test_value(model6.eval(dir, energy, time), 0.0, "Test sigma() setter");
  2076 + model6.sigma(spectrum);
  2077 + test_assert(model6.sigma() != NULL, "Test sigma() setter");
  2078 + test_value(model6.eval(dir, energy, time), 1.0, "Test sigma() setter");
  2079 +
  2080 + // Return
  2081 + return;
  2082 +}
  2083 +
  2084 +
  2085 +/***********************************************************************//**
2032 2086 * @brief Checks handling of CTA event bin
2033 2087 ***************************************************************************/
2034 2088 void TestGCTAObservation::test_event_bin(void)
... ...
inst/cta/test/test_CTA.hpp
  View file @c4b0b2a
... ... @@ -125,6 +125,7 @@ public:
125 125 void test_model_cube_bgd(void);
126 126 void test_model_irf_bgd(void);
127 127 void test_model_aeff_bgd(void);
  128 + void test_spatial_gauss_spectrum(void);
128 129 };
129 130  
130 131  
... ...