Eirik Hatlen / gammalib

Commit ce238e9657dc2d5e59af5a3af284557b11268759

Authored by Jürgen Knödlseder
1 parent da41e063

Add GCTAModelSpatialLookup class to handle CTA background lookups (#2878)

Inline Side-by-side
Showing 13 changed files with 1472 additions and 3 deletions
ChangeLog
  View file @ce238e9
1   -2019-04-13
  1 +2019-04-26
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
5 5  
  6 + Add GCTAModelSpatialLookup class to handle CTA background lookups (#2878)
6 7 Implement energy-dependent Gaussian CTA spatial background model (#2868)
7 8 Take instrument scale factors in On/Off analysis into account (#2867)
8 9 GWcs::wcs_set() method now sets radesys and equinox members (#2866)
... ...
1 1 New Features and Important Changes in GammaLib 1.6.0
2 2  
3   -13 April 2019
  3 +26 April 2019
4 4  
5 5  
6 6 1. Introduction
... ... @@ -21,6 +21,7 @@ The following classes have been added:
21 21 - GCTAModelSpatialMultiplicative
22 22 - GCTAModelSpatialRegistry
23 23 - GCTAModelSpatialGaussSpectrum
  24 +- GCTAModelSpatialLookup
24 25  
25 26 The following classes have been removed:
26 27 - GCTACubeSourceDiffuse
... ... @@ -332,6 +333,9 @@ unit tests for that response functions were added.
332 333  
333 334 20. CTA interface
334 335 -----------------
  336 +Add the GCTAModelSpatialLookup class that implements a lookup table for the
  337 +energy-dependent spatial modelling of the background (#2878).
  338 +
335 339 Add the GCTAModelSpatialGaussSpectrum class that implements an energy-dependent
336 340 Gaussian spatial background model (#2868).
337 341  
... ...
doc/source/admin/release_history/1.6.rst
  View file @ce238e9
... ... @@ -55,6 +55,8 @@ Bug fixes
55 55 Improvements
56 56 ------------
57 57  
  58 +* [`2878 <https://cta-redmine.irap.omp.eu/issues/2878>`_] -
  59 + Add ``GCTAModelSpatialLookup`` class to handle CTA background lookups
58 60 * [`2868 <https://cta-redmine.irap.omp.eu/issues/2868>`_] -
59 61 Implement energy-dependent Gaussian CTA spatial background model
60 62 * [`2867 <https://cta-redmine.irap.omp.eu/issues/2867>`_] -
... ...
inst/cta/Makefile.am
  View file @ce238e9
... ... @@ -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/GCTAModelSpatialLookup.cpp \
78 79 src/GCTAModelSpatialGaussSpectrum.cpp \
79 80 src/GCTAModelSpatialGradient.cpp \
80 81 src/GCTAModelSpatialMultiplicative.cpp \
... ... @@ -130,6 +131,7 @@ pkginclude_HEADERS = include/GCTATypemaps.hpp \
130 131 include/GCTACubeSource.hpp \
131 132 include/GCTAModelBackground.hpp \
132 133 include/GCTAModelSpatial.hpp \
  134 + include/GCTAModelSpatialLookup.hpp \
133 135 include/GCTAModelSpatialGaussSpectrum.hpp \
134 136 include/GCTAModelSpatialGradient.hpp \
135 137 include/GCTAModelSpatialMultiplicative.hpp \
... ...
inst/cta/include/GCTALib.hpp
  View file @ce238e9
... ... @@ -73,6 +73,7 @@
73 73 #include "GCTACubeSource.hpp"
74 74 #include "GCTAModelBackground.hpp"
75 75 #include "GCTAModelSpatial.hpp"
  76 +#include "GCTAModelSpatialLookup.hpp"
76 77 #include "GCTAModelSpatialGaussSpectrum.hpp"
77 78 #include "GCTAModelSpatialGradient.hpp"
78 79 #include "GCTAModelSpatialMultiplicative.hpp"
... ...
inst/cta/include/GCTAModelSpatialLookup.hpp 0 → 100644
  View file @ce238e9
  1 +/***************************************************************************
  2 + * GCTAModelSpatialLookup.hpp - Spatial lookup table model *
  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 GCTAModelSpatialLookup.hpp
  23 + * @brief Spatial lookup table model interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCTAMODELSPATIALLOOKUP_HPP
  28 +#define GCTAMODELSPATIALLOOKUP_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GCTAModelSpatial.hpp"
  33 +#include "GCTAResponseTable.hpp"
  34 +
  35 +/* __ Forward declarations _______________________________________________ */
  36 +class GFilename;
  37 +class GXmlElement;
  38 +class GCTAInstDir;
  39 +class GEnergy;
  40 +class GTime;
  41 +class GEbounds;
  42 +class GFitsTable;
  43 +class GFitsBinTable;
  44 +class GObservations;
  45 +class GCTAObservation;
  46 +
  47 +/* __ Constants __________________________________________________________ */
  48 +namespace gammalib {
  49 + const std::string extname_cta_spatial_lookup = "RADIAL BACKGROUND LOOKUP";
  50 +}
  51 +
  52 +
  53 +/***********************************************************************//**
  54 + * @class GCTAModelSpatialLookup
  55 + *
  56 + * @brief Spatial lookup table model class
  57 + ***************************************************************************/
  58 +class GCTAModelSpatialLookup : public GCTAModelSpatial {
  59 +
  60 +public:
  61 + // Constructors and destructors
  62 + GCTAModelSpatialLookup(void);
  63 + explicit GCTAModelSpatialLookup(const GFilename& filename);
  64 + explicit GCTAModelSpatialLookup(const GCTAResponseTable& table);
  65 + explicit GCTAModelSpatialLookup(const GXmlElement& xml);
  66 + GCTAModelSpatialLookup(const double& maxrad,
  67 + const double& radbin,
  68 + const GEbounds& ebds);
  69 + GCTAModelSpatialLookup(const GCTAModelSpatialLookup& model);
  70 + virtual ~GCTAModelSpatialLookup(void);
  71 +
  72 + // Operators
  73 + virtual GCTAModelSpatialLookup& operator=(const GCTAModelSpatialLookup& model);
  74 +
  75 + // Implemented pure virtual methods
  76 + virtual void clear(void);
  77 + virtual GCTAModelSpatialLookup* clone(void) const;
  78 + virtual std::string classname(void) const;
  79 + virtual std::string type(void) const;
  80 + virtual double eval(const GCTAInstDir& dir,
  81 + const GEnergy& energy,
  82 + const GTime& time,
  83 + const bool& gradients = false) const;
  84 + virtual double mc_max_value(const GCTAObservation& obs) const;
  85 + virtual void read(const GXmlElement& xml);
  86 + virtual void write(GXmlElement& xml) const;
  87 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  88 +
  89 + // Other methods
  90 + void fill(const GCTAObservation& obs);
  91 + void fill(const GObservations& obs);
  92 + const GCTAResponseTable& table(void) const;
  93 + void table(const GCTAResponseTable& table);
  94 + void read(const GFitsTable& table);
  95 + void write(GFitsBinTable& table) const;
  96 + void load(const GFilename& filename);
  97 + void save(const GFilename& filename,
  98 + const bool& clobber = false) const;
  99 + double norm(void) const;
  100 + void norm(const double& norm);
  101 +
  102 +protected:
  103 + // Protected methods
  104 + void init_members(void);
  105 + void copy_members(const GCTAModelSpatialLookup& model);
  106 + void free_members(void);
  107 + void prepare_table(void);
  108 + void normalise_table(void);
  109 + int table_index(const int& ienergy, const int& itheta) const;
  110 + void fill_buffer(const GCTAObservation& obs, std::vector<double>& buffer);
  111 + void set_from_buffer(const std::vector<double>& buffer);
  112 +
  113 + // Protected members
  114 + GFilename m_filename; //!< Name of lookup table
  115 + GCTAResponseTable m_lookup; //!< Lookup table
  116 + GModelPar m_norm; //!< Normalization factor
  117 + int m_inx_energy; //!< Energy index
  118 + int m_inx_theta; //!< Theta index
  119 +};
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Return class name
  124 + *
  125 + * @return String containing the class name ("GCTAModelSpatialLookup").
  126 + ***************************************************************************/
  127 +inline
  128 +std::string GCTAModelSpatialLookup::classname(void) const
  129 +{
  130 + return ("GCTAModelSpatialLookup");
  131 +}
  132 +
  133 +
  134 +/***********************************************************************//**
  135 + * @brief Return model type
  136 + *
  137 + * @return Model type "LookupTable".
  138 + ***************************************************************************/
  139 +inline
  140 +std::string GCTAModelSpatialLookup::type(void) const
  141 +{
  142 + return ("LookupTable");
  143 +}
  144 +
  145 +
  146 +/***********************************************************************//**
  147 + * @brief Return maximum function value for Monte Carlo simulations
  148 + *
  149 + * @param[in] obs CTA Observation.
  150 + * @return Maximum function value for Monte Carlo simulations.
  151 + *
  152 + * This method always returns 1.
  153 + ***************************************************************************/
  154 +inline
  155 +double GCTAModelSpatialLookup::mc_max_value(const GCTAObservation& obs) const
  156 +{
  157 + return 1.0;
  158 +}
  159 +
  160 +
  161 +/***********************************************************************//**
  162 + * @brief Return lookup table
  163 + *
  164 + * @return Lookup table.
  165 + *
  166 + * Returns the lookup table.
  167 + ***************************************************************************/
  168 +inline
  169 +const GCTAResponseTable& GCTAModelSpatialLookup::table(void) const
  170 +{
  171 + return m_lookup;
  172 +}
  173 +
  174 +
  175 +/***********************************************************************//**
  176 + * @brief Get lookup table model normalisation
  177 + *
  178 + * @return Lookup table model normalisation.
  179 + *
  180 + * Returns the normalisation of the lookup table model.
  181 + ***************************************************************************/
  182 +inline
  183 +double GCTAModelSpatialLookup::norm(void) const
  184 +{
  185 + return (m_norm.value());
  186 +}
  187 +
  188 +
  189 +/***********************************************************************//**
  190 + * @brief Set lookup table model normalisation
  191 + *
  192 + * @param[in] norm Lookup table model normalisation.
  193 + *
  194 + * Set the normalisation of the lookup table model.
  195 + ***************************************************************************/
  196 +inline
  197 +void GCTAModelSpatialLookup::norm(const double& norm)
  198 +{
  199 + m_norm.value(norm);
  200 + return;
  201 +}
  202 +
  203 +#endif /* GCTAMODELSPATIALLOOKUP_HPP */
... ...
inst/cta/pyext/GCTAModelSpatialLookup.i 0 → 100644
  View file @ce238e9
  1 +/***************************************************************************
  2 + * GCTAModelSpatialLookup.i - Spatial lookup table model *
  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 GCTAModelSpatialLookup.i
  23 + * @brief Spatial lookup table model interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCTAModelSpatialLookup.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCTAModelSpatialLookup
  34 + *
  35 + * @brief Spatial lookup table model class
  36 + ***************************************************************************/
  37 +class GCTAModelSpatialLookup : public GCTAModelSpatial {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + GCTAModelSpatialLookup(void);
  42 + explicit GCTAModelSpatialLookup(const GFilename& filename);
  43 + explicit GCTAModelSpatialLookup(const GCTAResponseTable& table);
  44 + explicit GCTAModelSpatialLookup(const GXmlElement& xml);
  45 + GCTAModelSpatialLookup(const double& radmax,
  46 + const double& radbin,
  47 + const GEbounds& ebds);
  48 + GCTAModelSpatialLookup(const GCTAModelSpatialLookup& model);
  49 + virtual ~GCTAModelSpatialLookup(void);
  50 +
  51 + // Operators
  52 + virtual GCTAModelSpatialLookup& operator=(const GCTAModelSpatialLookup& model);
  53 +
  54 + // Implemented pure virtual methods
  55 + virtual void clear(void);
  56 + virtual GCTAModelSpatialLookup* clone(void) const;
  57 + virtual std::string classname(void) const;
  58 + virtual std::string type(void) const;
  59 + virtual double eval(const GCTAInstDir& dir,
  60 + const GEnergy& energy,
  61 + const GTime& time,
  62 + const bool& gradients = false) const;
  63 + virtual double mc_max_value(const GCTAObservation& obs) const;
  64 + virtual void read(const GXmlElement& xml);
  65 + virtual void write(GXmlElement& xml) const;
  66 +
  67 + // Other methods
  68 + void fill(const GCTAObservation& obs);
  69 + void fill(const GObservations& obs);
  70 + const GCTAResponseTable& table(void) const;
  71 + void table(const GCTAResponseTable& table);
  72 + void read(const GFitsTable& table);
  73 + void write(GFitsBinTable& table) const;
  74 + void load(const GFilename& filename);
  75 + void save(const GFilename& filename,
  76 + const bool& clobber = false) const;
  77 + double norm(void) const;
  78 + void norm(const double& norm);
  79 +};
  80 +
  81 +
  82 +/***********************************************************************//**
  83 + * @brief GCTAModelSpatialLookup class extension
  84 + ***************************************************************************/
  85 +%extend GCTAModelSpatialLookup {
  86 + GCTAModelSpatialLookup copy() {
  87 + return (*self);
  88 + }
  89 +%pythoncode {
  90 + def __getstate__(self):
  91 + xml = gammalib.GXmlElement()
  92 + self.write(xml)
  93 + state = (xml,)
  94 + return state
  95 + def __setstate__(self, state):
  96 + if state[0].elements('parameter') == 0:
  97 + self.__init__()
  98 + else:
  99 + self.__init__(state[0])
  100 +}
  101 +};
... ...
inst/cta/pyext/cta.i
  View file @ce238e9
... ... @@ -308,6 +308,7 @@
308 308 %include "GCTAModelIrfBackground.i"
309 309 %include "GCTAModelAeffBackground.i"
310 310 %include "GCTAModelSpatial.i"
  311 +%include "GCTAModelSpatialLookup.i"
311 312 %include "GCTAModelSpatialGaussSpectrum.i"
312 313 %include "GCTAModelSpatialGradient.i"
313 314 %include "GCTAModelSpatialMultiplicative.i"
... ...
inst/cta/src/GCTAModelSpatialLookup.cpp 0 → 100644
  View file @ce238e9
  1 +/***************************************************************************
  2 + * GCTAModelSpatialLookup.cpp - Spatial lookup table model *
  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 GCTAModelSpatialLookup.cpp
  23 + * @brief Spatial lookup table model 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 "GEnergy.hpp"
  35 +#include "GTime.hpp"
  36 +#include "GEbounds.hpp"
  37 +#include "GFits.hpp"
  38 +#include "GFitsTable.hpp"
  39 +#include "GFitsBinTable.hpp"
  40 +#include "GObservations.hpp"
  41 +#include "GCTAObservation.hpp"
  42 +#include "GCTAInstDir.hpp"
  43 +#include "GCTAModelSpatialLookup.hpp"
  44 +#include "GCTAModelSpatialRegistry.hpp"
  45 +
  46 +/* __ Constants __________________________________________________________ */
  47 +
  48 +/* __ Globals ____________________________________________________________ */
  49 +const GCTAModelSpatialLookup g_cta_spatial_lookup_seed;
  50 +const GCTAModelSpatialRegistry g_cta_spatial_lookup_registry(&g_cta_spatial_lookup_seed);
  51 +
  52 +/* __ Method name definitions ____________________________________________ */
  53 +#define G_CONSTRUCTOR "GCTAModelSpatialLookup(double&, double&, GEbounds&)"
  54 +#define G_READ_XML "GCTAModelSpatialLookup::read(GXmlElement&)"
  55 +#define G_WRITE_XML "GCTAModelSpatialLookup::write(GXmlElement&)"
  56 +#define G_FILL_1 "GCTAModelSpatialLookup::fill(GCTAObservation&)"
  57 +#define G_FILL_2 "GCTAModelSpatialLookup::fill(GObservations&)"
  58 +#define G_READ_TABLE "GCTAModelSpatialLookup::read(GFitsTable&)"
  59 +#define G_LOAD "GCTAModelSpatialLookup::load(GFilename&)"
  60 +#define G_FILL_BUFFER "GCTAModelSpatialLookup::fill_buffer(GCTAObservation&,"\
  61 + " std::vector<double>&)"
  62 +
  63 +/* __ Macros _____________________________________________________________ */
  64 +
  65 +/* __ Coding definitions _________________________________________________ */
  66 +
  67 +/* __ Debug definitions __________________________________________________ */
  68 +
  69 +
  70 +/*==========================================================================
  71 + = =
  72 + = Constructors/destructors =
  73 + = =
  74 + ==========================================================================*/
  75 +
  76 +/***********************************************************************//**
  77 + * @brief Void constructor
  78 + ***************************************************************************/
  79 +GCTAModelSpatialLookup::GCTAModelSpatialLookup(void) : GCTAModelSpatial()
  80 +{
  81 + // Initialise members
  82 + init_members();
  83 +
  84 + // Return
  85 + return;
  86 +}
  87 +
  88 +
  89 +/***********************************************************************//**
  90 + * @brief Lookup table file constructor
  91 + *
  92 + * @param[in] filename Lookup table file name.
  93 + *
  94 + * Creates instance of a spatial lookup table model from a lookup table FITS
  95 + * file.
  96 + ***************************************************************************/
  97 +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const GFilename& filename) :
  98 + GCTAModelSpatial()
  99 +{
  100 + // Initialise members
  101 + init_members();
  102 +
  103 + // Load lookup table
  104 + load(filename);
  105 +
  106 + // Return
  107 + return;
  108 +}
  109 +
  110 +
  111 +/***********************************************************************//**
  112 + * @brief Response table constructor
  113 + *
  114 + * @param[in] table Response table.
  115 + *
  116 + * Creates instance of a spatial lookup table model from a response table.
  117 + ***************************************************************************/
  118 +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const GCTAResponseTable& table) :
  119 + GCTAModelSpatial()
  120 +{
  121 + // Initialise members
  122 + init_members();
  123 +
  124 + // Set lookup table
  125 + this->table(table);
  126 +
  127 + // Return
  128 + return;
  129 +}
  130 +
  131 +
  132 +/***********************************************************************//**
  133 + * @brief XML constructor
  134 + *
  135 + * @param[in] xml XML element.
  136 + *
  137 + * Creates instance of a spatial lookup table model from an XML element. See
  138 + * the read() method for information about the expected structure of the XML
  139 + * element.
  140 + ***************************************************************************/
  141 +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const GXmlElement& xml) :
  142 + GCTAModelSpatial()
  143 +{
  144 + // Initialise members
  145 + init_members();
  146 +
  147 + // Read information from XML element
  148 + read(xml);
  149 +
  150 + // Return
  151 + return;
  152 +}
  153 +
  154 +
  155 +/***********************************************************************//**
  156 + * @brief Lookup table constructor
  157 + *
  158 + * @param[in] maxrad Maximum radius (deg).
  159 + * @param[in] radbin Radial binsize (deg).
  160 + * @param[in] ebds Energy boundaries.
  161 + *
  162 + * @exception GException::invalid_argument
  163 + * Non positive radial binsize specified.
  164 + *
  165 + * Creates empty instance of a spatial lookup table model using the maximum
  166 + * lookup radius @p maxrad, the radial lookup binsize @p radbin, and some
  167 + * energy boundaries @p ebds. Following construction, the instance may be
  168 + * filled using the fill() methods.
  169 + *
  170 + * The method throws an exception if a non-positive radial bin size is
  171 + * specified.
  172 + ***************************************************************************/
  173 +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const double& maxrad,
  174 + const double& radbin,
  175 + const GEbounds& ebds) :
  176 + GCTAModelSpatial()
  177 +{
  178 + // Throw an exception if the radian binsize is not positive
  179 + if (radbin <= 0.0) {
  180 + std::string msg = "Non-positive radial bin size "+
  181 + gammalib::str(radbin)+" specified.";
  182 + throw GException::invalid_argument(G_CONSTRUCTOR, msg);
  183 + }
  184 +
  185 + // Initialise members
  186 + init_members();
  187 +
  188 + // Setup energy axis vectors
  189 + std::vector<double> eng_lo(ebds.size());
  190 + std::vector<double> eng_hi(ebds.size());
  191 + for (int i = 0; i < ebds.size(); ++i) {
  192 + eng_lo[i] = ebds.emin(i).TeV();
  193 + eng_hi[i] = ebds.emax(i).TeV();
  194 + }
  195 +
  196 + // Setup radial axis vector
  197 + int nrad = maxrad/radbin;
  198 + std::vector<double> rad_lo(nrad);
  199 + std::vector<double> rad_hi(nrad);
  200 + for (int i = 0; i < nrad; ++i) {
  201 + rad_lo[i] = double(i)*radbin;
  202 + rad_hi[i] = double(i+1)*radbin;
  203 + }
  204 +
  205 + // Setup response table
  206 + m_lookup.append_axis(eng_lo, eng_hi, "ENERG", "TeV");
  207 + m_lookup.append_axis(rad_lo, rad_hi, "THETA", "deg");
  208 + m_lookup.append_table("BKG", "1/sr");
  209 +
  210 + // Prepare table
  211 + prepare_table();
  212 +
  213 + // Return
  214 + return;
  215 +}
  216 +
  217 +
  218 +/***********************************************************************//**
  219 + * @brief Copy constructor
  220 + *
  221 + * @param[in] model Lookup table spatial model.
  222 + ***************************************************************************/
  223 +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const GCTAModelSpatialLookup& model) :
  224 + GCTAModelSpatial(model)
  225 +{
  226 + // Initialise members
  227 + init_members();
  228 +
  229 + // Copy members
  230 + copy_members(model);
  231 +
  232 + // Return
  233 + return;
  234 +}
  235 +
  236 +
  237 +/***********************************************************************//**
  238 + * @brief Destructor
  239 + ***************************************************************************/
  240 +GCTAModelSpatialLookup::~GCTAModelSpatialLookup(void)
  241 +{
  242 + // Free members
  243 + free_members();
  244 +
  245 + // Return
  246 + return;
  247 +}
  248 +
  249 +
  250 +/*==========================================================================
  251 + = =
  252 + = Operators =
  253 + = =
  254 + ==========================================================================*/
  255 +
  256 +/***********************************************************************//**
  257 + * @brief Assignment operator
  258 + *
  259 + * @param[in] model Lookup table spatial model.
  260 + ***************************************************************************/
  261 +GCTAModelSpatialLookup& GCTAModelSpatialLookup::operator=(const GCTAModelSpatialLookup& model)
  262 +{
  263 + // Execute only if object is not identical
  264 + if (this != &model) {
  265 +
  266 + // Copy base class members
  267 + this->GCTAModelSpatial::operator=(model);
  268 +
  269 + // Free members
  270 + free_members();
  271 +
  272 + // Initialise members
  273 + init_members();
  274 +
  275 + // Copy members
  276 + copy_members(model);
  277 +
  278 + } // endif: object was not identical
  279 +
  280 + // Return
  281 + return *this;
  282 +}
  283 +
  284 +
  285 +/*==========================================================================
  286 + = =
  287 + = Public methods =
  288 + = =
  289 + ==========================================================================*/
  290 +
  291 +/***********************************************************************//**
  292 + * @brief Clear instance
  293 + ***************************************************************************/
  294 +void GCTAModelSpatialLookup::clear(void)
  295 +{
  296 + // Free class members (base and derived classes, derived class first)
  297 + free_members();
  298 + this->GCTAModelSpatial::free_members();
  299 +
  300 + // Initialise members
  301 + this->GCTAModelSpatial::init_members();
  302 + init_members();
  303 +
  304 + // Return
  305 + return;
  306 +}
  307 +
  308 +
  309 +/***********************************************************************//**
  310 + * @brief Clone instance
  311 + ***************************************************************************/
  312 +GCTAModelSpatialLookup* GCTAModelSpatialLookup::clone(void) const
  313 +{
  314 + return new GCTAModelSpatialLookup(*this);
  315 +}
  316 +
  317 +
  318 +/***********************************************************************//**
  319 + * @brief Evaluate function
  320 + *
  321 + * @param[in] dir Event direction.
  322 + * @param[in] energy Event energy.
  323 + * @param[in] time Event time.
  324 + * @param[in] gradients Compute gradients?
  325 + * @return Function value
  326 + *
  327 + * Evaluates the lookup table model for a given event direction and energy.
  328 + * The evaluation is done using bi-linear interpolation in the logarithm
  329 + * of the energy and the offset angle radius. The lookup table value is
  330 + * multiplied with a normalisation parameter. The method always returns a
  331 + * non-negative result.
  332 + *
  333 + * If the @p gradients flag is true the method will also compute the partial
  334 + * derivatives of the normalisation parameter.
  335 + ***************************************************************************/
  336 +double GCTAModelSpatialLookup::eval(const GCTAInstDir& dir,
  337 + const GEnergy& energy,
  338 + const GTime& time,
  339 + const bool& gradients) const
  340 +{
  341 + // Compute offset angle in radians
  342 + double offset = dir.theta();
  343 +
  344 + // Determine function value by bilinear interpolation
  345 + double value = ((m_lookup.elements() > 0) &&
  346 + (m_lookup.tables() > 0))
  347 + ? m_lookup(0, energy.log10TeV(), offset) : 0.0;
  348 +
  349 + // Make sure that function does not become negative
  350 + if (value < 0.0) {
  351 + value = 0.0;
  352 + }
  353 +
  354 + // Optionally compute partial derivatives
  355 + if (gradients) {
  356 +
  357 + // Compute partial derivative of the normalisation parameter
  358 + double g_norm = (m_norm.is_free()) ? value * m_norm.scale() : 0.0;
  359 +
  360 + // Set gradient
  361 + m_norm.factor_gradient(g_norm);
  362 +
  363 + } // endif: computed partial derivatives
  364 +
  365 + // Return intensity times normalization factor
  366 + return (value * m_norm.value());
  367 +}
  368 +
  369 +
  370 +/***********************************************************************//**
  371 + * @brief Read model from XML element
  372 + *
  373 + * @param[in] xml XML element.
  374 + *
  375 + * Read the lookup table model information from an XML element. The XML
  376 + * element needs to be of the following format:
  377 + *
  378 + * <spatialModel type="LookupTable" file="lookuptable.fits">
  379 + * <parameter name="Normalization" ../>
  380 + * </spatialModel>
  381 + *
  382 + * The @p file attribute provides the filename of the lookup table FITS file.
  383 + * The filename may be either an absolute filename (starting with '/') or a
  384 + * relative filename. If no access path is given, the file is expected to
  385 + * reside in the same location as the XML file.
  386 + ***************************************************************************/
  387 +void GCTAModelSpatialLookup::read(const GXmlElement& xml)
  388 +{
  389 + // Clear model
  390 + clear();
  391 +
  392 + // Load lookup table
  393 + load(gammalib::xml_file_expand(xml, xml.attribute("file")));
  394 +
  395 + // Get parameter pointers
  396 + const GXmlElement* norm = gammalib::xml_get_par(G_READ_XML, xml, m_norm.name());
  397 +
  398 + // Read parameters
  399 + m_norm.read(*norm);
  400 +
  401 + // Return
  402 + return;
  403 +}
  404 +
  405 +
  406 +/***********************************************************************//**
  407 + * @brief Write model into XML element
  408 + *
  409 + * @param[in] xml XML element.
  410 + *
  411 + * @exception GException::invalid_value
  412 + * Spatial model is not of valid type.
  413 + *
  414 + * Write the lookup table model information into an XML element. The XML
  415 + * element will be of the following format:
  416 + *
  417 + * <spatialModel type="LookupTable" file="lookuptable.fits">
  418 + * <parameter name="Normalization" ../>
  419 + * </spatialModel>
  420 + *
  421 + * The @p file attribute provides the filename of the lookup table FITS file.
  422 + * The filename may be either an absolute filename (starting with '/') or a
  423 + * relative filename. If no access path is given, the file is expected to
  424 + * reside in the same location as the XML file.
  425 + ***************************************************************************/
  426 +void GCTAModelSpatialLookup::write(GXmlElement& xml) const
  427 +{
  428 + // Set model type
  429 + if (xml.attribute("type") == "") {
  430 + xml.attribute("type", type());
  431 + }
  432 +
  433 + // Verify model type
  434 + if (xml.attribute("type") != type()) {
  435 + std::string msg = "Spatial model \""+xml.attribute("type")+
  436 + "\" is not of type \""+type()+"\".";
  437 + throw GException::invalid_value(G_WRITE_XML, msg);
  438 + }
  439 +
  440 + // Set lookup table file name
  441 + xml.attribute("file", gammalib::xml_file_reduce(xml, m_filename));
  442 +
  443 + // Get XML parameters
  444 + GXmlElement* norm = gammalib::xml_need_par(G_WRITE_XML, xml, m_norm.name());
  445 +
  446 + // Write parameters
  447 + m_norm.write(*norm);
  448 +
  449 + // Return
  450 + return;
  451 +}
  452 +
  453 +
  454 +/***********************************************************************//**
  455 + * @brief Fill lookup table with events from one CTA observation
  456 + *
  457 + * @param[in] obs CTA observation.
  458 + *
  459 + * @exception GException::invalid_value
  460 + * No lookup table has been defined.
  461 + *
  462 + * Fills lookup table with events from one CTA observation. The method needs
  463 + * that a lookup table is defined before. In case that no lookup table is
  464 + * defined the method throws an exception.
  465 + ***************************************************************************/
  466 +void GCTAModelSpatialLookup::fill(const GCTAObservation& obs)
  467 +{
  468 + // Throw an exception if response table is not initialised
  469 + if (m_lookup.elements() == 0) {
  470 + std::string msg = "Attempting to fill a yet undefined lookup table. "
  471 + "Please allocate a valid lookup table before "
  472 + "calling this method.";
  473 + throw GException::invalid_value(G_FILL_1, msg);
  474 + }
  475 + if (m_lookup.tables() == 0) {
  476 + std::string msg = "Attempting to fill non-existing lookup table. "
  477 + "While the lookup table axes were defined, no table "
  478 + "exists. Please allocate a valid lookup table before "
  479 + "calling this method.";
  480 + throw GException::invalid_value(G_FILL_1, msg);
  481 + }
  482 +
  483 + // Allocate buffer for events
  484 + std::vector<double> counts(m_lookup.elements());
  485 +
  486 + // Fill counts into buffer
  487 + fill_buffer(obs, counts);
  488 +
  489 + // Set lookup table from buffer
  490 + set_from_buffer(counts);
  491 +
  492 + // Return
  493 + return;
  494 +}
  495 +
  496 +
  497 +/***********************************************************************//**
  498 + * @brief Fill lookup table with CTA events in observation container
  499 + *
  500 + * @param[in] obs Observation container.
  501 + *
  502 + * @exception GException::invalid_value
  503 + * No lookup table has been defined.
  504 + *
  505 + * Fills lookup table with CTA events in observation container. The method
  506 + * needs that a lookup table is defined before. In case that no lookup table
  507 + * is defined the method throws an exception.
  508 + ***************************************************************************/
  509 +void GCTAModelSpatialLookup::fill(const GObservations& obs)
  510 +{
  511 + // Throw an exception if response table is not initialised
  512 + if (m_lookup.elements() == 0) {
  513 + std::string msg = "Attempting to fill a yet undefined lookup table. "
  514 + "Please allocate a valid lookup table before "
  515 + "calling this method.";
  516 + throw GException::invalid_value(G_FILL_2, msg);
  517 + }
  518 + if (m_lookup.tables() == 0) {
  519 + std::string msg = "Attempting to fill non-existing lookup table. "
  520 + "While the lookup table axes were defined, no table "
  521 + "exists. Please allocate a valid lookup table before "
  522 + "calling this method.";
  523 + throw GException::invalid_value(G_FILL_2, msg);
  524 + }
  525 +
  526 + // Allocate buffer for events
  527 + std::vector<double> counts(m_lookup.elements());
  528 +
  529 + // Loop over all observations in container
  530 + for (int i = 0; i < obs.size(); ++i) {
  531 +
  532 + // Get observation and continue only if it is a CTA observation
  533 + const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
  534 +
  535 + // Skip observation if it's not CTA
  536 + if (cta == NULL) {
  537 + continue;
  538 + }
  539 +
  540 + // Skip observation if we have a binned observation
  541 + if (cta->eventtype() == "CountsCube") {
  542 + continue;
  543 + }
  544 +
  545 + // Fill counts into buffer
  546 + fill_buffer(*cta, counts);
  547 +
  548 + }
  549 +
  550 + // Set lookup table from buffer
  551 + set_from_buffer(counts);
  552 +
  553 + // Return
  554 + return;
  555 +}
  556 +
  557 +
  558 +/***********************************************************************//**
  559 + * @brief Assign lookup table
  560 + *
  561 + * @param[in] table Lookup table.
  562 + *
  563 + * Assigns lookup table.
  564 + ***************************************************************************/
  565 +void GCTAModelSpatialLookup::table(const GCTAResponseTable& table)
  566 +{
  567 + // Assign lookup table
  568 + m_lookup = table;
  569 +
  570 + // Prepare table
  571 + prepare_table();
  572 +
  573 + // Return
  574 + return;
  575 +}
  576 +
  577 +
  578 +/***********************************************************************//**
  579 + * @brief Read lookup table from FITS table
  580 + *
  581 + * @param[in] table FITS table.
  582 + *
  583 + * @exception GException::invalid_value
  584 + * Response table is not two-dimensional.
  585 + *
  586 + * Reads the lookup table form the FITS @p table. The following column names
  587 + * are mandatory:
  588 + *
  589 + * ENERG_LO - Energy lower bin boundaries
  590 + * ENERG_HI - Energy upper bin boundaries
  591 + * THETA_LO - Offset angle lower bin boundaries
  592 + * THETA_HI - Offset angle upper bin boundaries
  593 + *
  594 + * After reading, the method assures that the lookup table is properly
  595 + * normalised, i.e. that for each energy vector the radial profile has a
  596 + * maximum value of one.
  597 + ***************************************************************************/
  598 +void GCTAModelSpatialLookup::read(const GFitsTable& table)
  599 +{
  600 + // Clear lookup table
  601 + m_lookup.clear();
  602 +
  603 + // Read lookup table
  604 + m_lookup.read(table);
  605 +
  606 + // Throw an exception if the table is not two-dimensional
  607 + if (m_lookup.axes() != 2) {
  608 + std::string msg = "Expected two-dimensional lookup table but found "+
  609 + gammalib::str(m_lookup.axes())+" dimensions. Please "
  610 + "specify a two-dimensional lookup table.";
  611 + throw GException::invalid_value(G_READ_TABLE, msg);
  612 + }
  613 +
  614 + // Prepare table
  615 + prepare_table();
  616 +
  617 + // Return
  618 + return;
  619 +}
  620 +
  621 +
  622 +/***********************************************************************//**
  623 + * @brief Write lookup table into FITS binary table
  624 + *
  625 + * @param[in] table FITS binary table.
  626 + *
  627 + * Writes lookup table into a FITS binary @p table.
  628 + ***************************************************************************/
  629 +void GCTAModelSpatialLookup::write(GFitsBinTable& table) const
  630 +{
  631 + // Write lookup table
  632 + m_lookup.write(table);
  633 +
  634 + // Return
  635 + return;
  636 +}
  637 +
  638 +
  639 +/***********************************************************************//**
  640 + * @brief Load lookup table
  641 + *
  642 + * @param[in] filename Load lookup table file name.
  643 + *
  644 + * Loads lookup table from FITS file.
  645 + ***************************************************************************/
  646 +void GCTAModelSpatialLookup::load(const GFilename& filename)
  647 +{
  648 + // Clear sigma spectrum
  649 + clear();
  650 +
  651 + // Continue only if filename is not empty
  652 + if (!filename.is_empty()) {
  653 +
  654 + // Open FITS file
  655 + GFits fits(filename);
  656 +
  657 + // Open FITS table
  658 + const GFitsTable* table = NULL;
  659 + if (filename.has_extname()) {
  660 + table = fits.table(filename.extname());
  661 + }
  662 + else if (filename.has_extno()) {
  663 + table = fits.table(filename.extno());
  664 + }
  665 + else {
  666 + table = fits.table(1);
  667 + }
  668 +
  669 + // If no FITS table is found then throw an exception
  670 + if (table == NULL) {
  671 + std::string msg = "No lookup table found in file \""+
  672 + filename.url()+"\".";
  673 + throw GException::file_error(G_LOAD, msg);
  674 + }
  675 +
  676 + // Read lookup table
  677 + read(*table);
  678 +
  679 + // Close FITS file
  680 + fits.close();
  681 +
  682 + // Store filename
  683 + m_filename = filename;
  684 +
  685 + } // endif: filename was not empty
  686 +
  687 + // Return
  688 + return;
  689 +}
  690 +
  691 +
  692 +/***********************************************************************//**
  693 + * @brief Save lookup table into FITS file
  694 + *
  695 + * @param[in] filename FITS file name.
  696 + * @param[in] clobber Overwrite existing file?
  697 + *
  698 + * Saves lookup table into a FITS file. If a file with the given @p filename
  699 + * does not yet exist it will be created, otherwise the method opens the
  700 + * existing file. The method will create a (or replace an existing)
  701 + * effective area extension. The extension name can be specified as part
  702 + * of the @p filename, or if no extension name is given, is assumed to be
  703 + * `RADIAL BACKGROUND LOOKUP`.
  704 + *
  705 + * An existing file will only be modified if the @p clobber flag is set to
  706 + * true.
  707 + ***************************************************************************/
  708 +void GCTAModelSpatialLookup::save(const GFilename& filename, const bool& clobber) const
  709 +{
  710 + // Get extension name
  711 + std::string extname = filename.extname(gammalib::extname_cta_spatial_lookup);
  712 +
  713 + // Open or create FITS file (without extension name since the requested
  714 + // extension may not yet exist in the file)
  715 + GFits fits(filename.url(), true);
  716 +
  717 + // Remove extension if it exists already
  718 + if (fits.contains(extname)) {
  719 + fits.remove(extname);
  720 + }
  721 +
  722 + // Create binary table
  723 + GFitsBinTable table;
  724 +
  725 + // Write the lookup table
  726 + write(table);
  727 +
  728 + // Set binary table extension name
  729 + table.extname(extname);
  730 +
  731 + // Append table to FITS file
  732 + fits.append(table);
  733 +
  734 + // Save to file
  735 + fits.save(clobber);
  736 +
  737 + // Return
  738 + return;
  739 +}
  740 +
  741 +
  742 +/***********************************************************************//**
  743 + * @brief Print lookup table information
  744 + *
  745 + * @param[in] chatter Chattiness.
  746 + * @return String containing lookup table information.
  747 + ***************************************************************************/
  748 +std::string GCTAModelSpatialLookup::print(const GChatter& chatter) const
  749 +{
  750 + // Initialise result string
  751 + std::string result;
  752 +
  753 + // Continue only if chatter is not silent
  754 + if (chatter != SILENT) {
  755 +
  756 + // Initialise information
  757 + int nebins = 0;
  758 + int nthetabins = 0;
  759 + double emin = 0.0;
  760 + double emax = 0.0;
  761 + double omin = 0.0;
  762 + double omax = 0.0;
  763 +
  764 + // Extract information if there are axes in the response table
  765 + if (m_lookup.axes() > 0) {
  766 + nebins = m_lookup.axis_bins(m_inx_energy);
  767 + nthetabins = m_lookup.axis_bins(m_inx_theta);
  768 + emin = m_lookup.axis_lo(m_inx_energy,0);
  769 + emax = m_lookup.axis_hi(m_inx_energy,nebins-1);
  770 + omin = m_lookup.axis_lo(m_inx_theta,0);
  771 + omax = m_lookup.axis_hi(m_inx_theta,nthetabins-1);
  772 + }
  773 +
  774 + // Append header
  775 + result.append("=== GCTAModelSpatialLookup ===");
  776 +
  777 + // Append information
  778 + result.append("\n"+gammalib::parformat("Filename")+m_filename);
  779 + result.append("\n"+gammalib::parformat("Number of energy bins") +
  780 + gammalib::str(nebins));
  781 + result.append("\n"+gammalib::parformat("Number of offset bins") +
  782 + gammalib::str(nthetabins));
  783 + result.append("\n"+gammalib::parformat("Energy range"));
  784 + result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
  785 + result.append("\n"+gammalib::parformat("Offset angle range"));
  786 + result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
  787 +
  788 + } // endif: chatter was not silent
  789 +
  790 + // Return result
  791 + return result;
  792 +}
  793 +
  794 +
  795 +/*==========================================================================
  796 + = =
  797 + = Private methods =
  798 + = =
  799 + ==========================================================================*/
  800 +
  801 +/***********************************************************************//**
  802 + * @brief Initialise class members
  803 + ***************************************************************************/
  804 +void GCTAModelSpatialLookup::init_members(void)
  805 +{
  806 + // Initialise members
  807 + m_filename.clear();
  808 + m_lookup.clear();
  809 + m_inx_energy = 0;
  810 + m_inx_theta = 1;
  811 +
  812 + // Initialise lookup table normalisation
  813 + m_norm.clear();
  814 + m_norm.name("Normalization");
  815 + m_norm.value(1.0);
  816 + m_norm.scale(1.0);
  817 + m_norm.range(0.001, 1000.0);
  818 + m_norm.gradient(0.0);
  819 + m_norm.fix();
  820 + m_norm.has_grad(true);
  821 +
  822 + // Set parameter pointer(s)
  823 + m_pars.clear();
  824 + m_pars.push_back(&m_norm);
  825 +
  826 + // Return
  827 + return;
  828 +}
  829 +
  830 +
  831 +/***********************************************************************//**
  832 + * @brief Copy class members
  833 + *
  834 + * @param[in] model Lookup table model.
  835 + ***************************************************************************/
  836 +void GCTAModelSpatialLookup::copy_members(const GCTAModelSpatialLookup& model)
  837 +{
  838 + // Copy members
  839 + m_filename = model.m_filename;
  840 + m_lookup = model.m_lookup;
  841 + m_norm = model.m_norm;
  842 + m_inx_energy = model.m_inx_energy;
  843 + m_inx_theta = model.m_inx_theta;
  844 +
  845 + // Set parameter pointer(s)
  846 + m_pars.clear();
  847 + m_pars.push_back(&m_norm);
  848 +
  849 + // Return
  850 + return;
  851 +}
  852 +
  853 +
  854 +/***********************************************************************//**
  855 + * @brief Delete class members
  856 + ***************************************************************************/
  857 +void GCTAModelSpatialLookup::free_members(void)
  858 +{
  859 + // Return
  860 + return;
  861 +}
  862 +
  863 +
  864 +/***********************************************************************//**
  865 + * @brief Prepare lookup table indices
  866 + *
  867 + * Prepares a lookup table for usage.
  868 + *
  869 + * The method sets the data members m_inx_energy and m_inx_theta which
  870 + * determine the axes order.
  871 + *
  872 + * The method furthermore sets the energy axis to logarithmic scale and the
  873 + * offset angles to radians.
  874 + *
  875 + * Finally, the method assures that the maximum values for each radial
  876 + * profile are one.
  877 + ***************************************************************************/
  878 +void GCTAModelSpatialLookup::prepare_table(void)
  879 +{
  880 + // Get mandatory indices (throw exception if not found)
  881 + m_inx_energy = m_lookup.axis("ENERG");
  882 + m_inx_theta = m_lookup.axis("THETA");
  883 +
  884 + // Set energy axis to logarithmic scale
  885 + m_lookup.axis_log10(m_inx_energy);
  886 +
  887 + // Set offset angle axis to radians
  888 + m_lookup.axis_radians(m_inx_theta);
  889 +
  890 + // Normalise radial profiles
  891 + normalise_table();
  892 +
  893 + // Return
  894 + return;
  895 +}
  896 +
  897 +
  898 +/***********************************************************************//**
  899 + * @brief Normalise lookup table
  900 + *
  901 + * Normalises lookup table so that each radial profile has a maximum value
  902 + * of one.
  903 + ***************************************************************************/
  904 +void GCTAModelSpatialLookup::normalise_table(void)
  905 +{
  906 + // Get axes dimensions
  907 + int energy_size = m_lookup.axis_bins(m_inx_energy);
  908 + int theta_size = m_lookup.axis_bins(m_inx_theta);
  909 +
  910 + // Loop over energies
  911 + for (int i_energy = 0; i_energy < energy_size; ++i_energy) {
  912 +
  913 + // Initialise maximum radial profile value
  914 + double rad_max = 0.0;
  915 +
  916 + // Determine maximum radial profile value
  917 + for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
  918 + int i = table_index(i_energy, i_theta);
  919 + if (m_lookup(0,i) > rad_max) {
  920 + rad_max = m_lookup(0,i);
  921 + }
  922 + }
  923 +
  924 + // If maximum value is positive then normalise profile
  925 + if (rad_max > 0.0) {
  926 + for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
  927 + int i = table_index(i_energy, i_theta);
  928 + m_lookup(0,i) /= rad_max;
  929 + }
  930 + }
  931 +
  932 + } // endfor: looped over all energies
  933 +
  934 + // Return
  935 + return;
  936 +}
  937 +
  938 +
  939 +/***********************************************************************//**
  940 + * @brief Return index of lookup table element
  941 + *
  942 + * @param[in] ienergy Energy index.
  943 + * @param[in] itheta Offset index.
  944 + * @return Index of lookup table element.
  945 + ***************************************************************************/
  946 +int GCTAModelSpatialLookup::table_index(const int& ienergy,
  947 + const int& itheta) const
  948 +{
  949 + // Set index vector
  950 + int inx[2];
  951 + inx[m_inx_energy] = ienergy;
  952 + inx[m_inx_theta] = itheta;
  953 +
  954 + // Compute index
  955 + int index = inx[0] + inx[1] * m_lookup.axis_bins(0);
  956 +
  957 + // Return index
  958 + return index;
  959 +}
  960 +
  961 +
  962 +/***********************************************************************//**
  963 + * @brief Fill buffer from events in CTA observation
  964 + *
  965 + * @param[in] obs CTA observation.
  966 + * @param[in,out] buffer Counts buffer.
  967 + ***************************************************************************/
  968 +void GCTAModelSpatialLookup::fill_buffer(const GCTAObservation& obs,
  969 + std::vector<double>& buffer)
  970 +{
  971 + // Make sure that the observation holds a CTA event list. If this is
  972 + // not the case then throw an exception.
  973 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  974 + if (events == NULL) {
  975 + std::string msg = "CTA Observation does not contain an event "
  976 + "list. An event list is needed to fill the "
  977 + "lookup table.";
  978 + throw GException::invalid_value(G_FILL_BUFFER, msg);
  979 + }
  980 +
  981 + // Get axes dimensions
  982 + int energy_size = m_lookup.axis_bins(m_inx_energy);
  983 + int theta_size = m_lookup.axis_bins(m_inx_theta);
  984 +
  985 + // Fill buffer
  986 + for (int i = 0; i < events->size(); ++i) {
  987 +
  988 + // Get event
  989 + const GCTAEventAtom* event = (*events)[i];
  990 +
  991 + // Determine offset angle in degrees
  992 + GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
  993 + double theta = inst->theta() * gammalib::rad2deg;
  994 +
  995 + // Determine offset angle bin
  996 + int i_theta = 0;
  997 + for (; i_theta < theta_size; ++i_theta) {
  998 + if ((theta >= m_lookup.axis_lo(m_inx_theta,i_theta)) &&
  999 + (theta <= m_lookup.axis_hi(m_inx_theta,i_theta))) {
  1000 + break;
  1001 + }
  1002 + }
  1003 +
  1004 + // Skip event if no valid offset angle bin was found
  1005 + if (i_theta >= theta_size) {
  1006 + continue;
  1007 + }
  1008 +
  1009 + // Determine event energy in TeV
  1010 + double energy = event->energy().TeV();
  1011 +
  1012 + // Determine energy bin
  1013 + int i_energy = 0;
  1014 + for (; i_energy < energy_size; ++i_energy) {
  1015 + if ((energy >= m_lookup.axis_lo(m_inx_energy,i_energy)) &&
  1016 + (energy <= m_lookup.axis_hi(m_inx_energy,i_energy))) {
  1017 + break;
  1018 + }
  1019 + }
  1020 +
  1021 + // Skip event if no valid offset angle bin was found
  1022 + if (i_energy >= theta_size) {
  1023 + continue;
  1024 + }
  1025 +
  1026 + // Fill event into bugger
  1027 + int inx = table_index(i_energy, i_theta);
  1028 + buffer[inx] += 1.0;
  1029 +
  1030 + } // endfor: looped over all events
  1031 +
  1032 + // Return
  1033 + return;
  1034 +}
  1035 +
  1036 +
  1037 +/***********************************************************************//**
  1038 + * @brief Set lookup table from buffer
  1039 + *
  1040 + * @param[in] buffer Counts buffer.
  1041 + ***************************************************************************/
  1042 +void GCTAModelSpatialLookup::set_from_buffer(const std::vector<double>& buffer)
  1043 +{
  1044 + // Get axes dimensions
  1045 + int energy_size = m_lookup.axis_bins(m_inx_energy);
  1046 + int theta_size = m_lookup.axis_bins(m_inx_theta);
  1047 +
  1048 + // Loop over offset angle
  1049 + for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
  1050 +
  1051 + // Compute solid angle for offset angle bin
  1052 + double theta_min = m_lookup.axis_lo(m_inx_theta,i_theta) * gammalib::deg2rad;
  1053 + double theta_max = m_lookup.axis_hi(m_inx_theta,i_theta) * gammalib::deg2rad;
  1054 + double area = gammalib::pi * (theta_max*theta_max - theta_min*theta_min);
  1055 +
  1056 + // Fill energy vector into lookup table
  1057 + for (int i_energy = 0; i_energy < energy_size; ++i_energy) {
  1058 + int inx = table_index(i_energy, i_theta);
  1059 + m_lookup(0,inx) = buffer[inx] / area;
  1060 + }
  1061 +
  1062 + } // endfor: looped over offset angles
  1063 +
  1064 + // Normalise radial profiles
  1065 + normalise_table();
  1066 +
  1067 + // Return
  1068 + return;
  1069 +}
... ...
inst/cta/test/data/bkg_lookup.fits 0 → 100644
View file @ce238e9
No preview for this file type
inst/cta/test/data/cta_model_bkg_lookup.xml 0 → 100644
  View file @ce238e9
  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="LookupTable" file="bkg_lookup.fits">
  10 + <parameter name="Normalization" scale="1" value="1" min="0.01" max="10.0" free="0"/>
  11 + </spatialModel>
  12 + </source>
  13 +</source_library>
... ...
inst/cta/test/test_CTA.cpp
  View file @ce238e9
... ... @@ -63,6 +63,8 @@ const std::string cta_bgd_gauss_e_xml = datadir+&quot;/cta_model_bgd_gauss_e.xml&quot;;
63 63 const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml";
64 64 const std::string cta_irf_bgd_xml = datadir+"/cta_model_irf_bgd.xml";
65 65 const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml";
  66 +const std::string cta_bkg_lookup = datadir+"/bkg_lookup.fits";
  67 +const std::string cta_bkg_lookup_xml = datadir+"/cta_model_bkg_lookup.xml";
66 68 const std::string cta_caldb_king = cta_caldb+"/data/cta/e/bcf/IFAE20120510_50h_King";
67 69 const std::string cta_irf_king = "irf_file.fits";
68 70 const std::string cta_psf_table = caldbdir+"/psf_table.fits[PSF_2D_TABLE]";
... ... @@ -70,7 +72,7 @@ const std::string cta_perf_table = caldbdir+&quot;/cta_dummy_irf.dat&quot;;
70 72 const std::string cta_edisp_rmf = caldbdir+"/dc1/rmf.fits";
71 73 const std::string cta_edisp_2D = caldbdir+"/edisp_matrix.fits";
72 74 const std::string cta_bgd_3D = caldbdir+"/edisp_matrix.fits";
73   -const std::string cta_modbck_fit = datadir+"/bg_test.fits";
  75 +//const std::string cta_modbck_fit = datadir+"/bg_test.fits";
74 76  
75 77 /* __ Test files for stacked analysis (based on Prod2::South_0.5h) _______ */
76 78 const std::string cta_stacked_xml = datadir+"/stacked_obs.xml";
... ... @@ -177,6 +179,8 @@ void TestGCTAModel::set(void)
177 179 "Test CTA IRF background model");
178 180 append(static_cast<pfunction>(&TestGCTAModel::test_model_aeff_bgd),
179 181 "Test CTA Aeff background model");
  182 + append(static_cast<pfunction>(&TestGCTAModel::test_spatial_lookup),
  183 + "Test GCTAModelSpatialLookup class");
180 184 append(static_cast<pfunction>(&TestGCTAModel::test_spatial_gauss_spectrum),
181 185 "Test GCTAModelSpatialGaussSpectrum class");
182 186  
... ... @@ -2032,6 +2036,73 @@ void TestGCTAModel::test_model_aeff_bgd(void)
2032 2036  
2033 2037  
2034 2038 /***********************************************************************//**
  2039 + * @brief Test GCTAModelSpatialLookup class
  2040 + ***************************************************************************/
  2041 +void TestGCTAModel::test_spatial_lookup(void)
  2042 +{
  2043 + // Set instrument direction, energy and time
  2044 + GCTAInstDir dir(0.0, 0.0);
  2045 + GEnergy energy(1.0, "TeV");
  2046 + GTime time;
  2047 +
  2048 + // Test void constructor
  2049 + GCTAModelSpatialLookup model1;
  2050 + test_assert(!model1.table().has_table("BKG"), "Test table() for void constructor");
  2051 + test_value(model1.table().elements(), 0, "Test table() for void constructor");
  2052 + test_value(model1.eval(dir, energy, time), 0.0, "Test eval() for void constructor");
  2053 + test_value(model1.norm(), 1.0, "Test norm() for void constructor");
  2054 +
  2055 + // Test filename constructor
  2056 + GCTAModelSpatialLookup model2(cta_bkg_lookup);
  2057 + test_assert(model2.table().has_table("BKG"), "Test table() for filename constructor");
  2058 + test_value(model2.table().elements(), 100, "Test table() for filename constructor");
  2059 + test_value(model2.eval(dir, energy, time), 1.0140199, "Test eval() for filename constructor");
  2060 + test_value(model2.norm(), 1.0, "Test norm() for filename constructor");
  2061 +
  2062 + // Test response table constructor
  2063 + GFits fits(cta_bkg_lookup);
  2064 + const GFitsTable& hdu = *fits.table(1);
  2065 + GCTAResponseTable table(hdu);
  2066 + GCTAModelSpatialLookup model3(table);
  2067 + test_assert(model3.table().has_table("BKG"), "Test table() for response table constructor");
  2068 + test_value(model3.table().elements(), 100, "Test table() for response table constructor");
  2069 + test_value(model3.eval(dir, energy, time), 1.0140199, "Test eval() for response table constructor");
  2070 + test_value(model3.norm(), 1.0, "Test norm() for response table constructor");
  2071 +
  2072 + // Test XML constructor
  2073 + GXml xml(cta_bkg_lookup_xml);
  2074 + GXmlElement* element = xml.element("source_library > source[0] > spatialModel");
  2075 + GCTAModelSpatialLookup model4(*element);
  2076 + test_assert(model4.table().has_table("BKG"), "Test table() for XML constructor");
  2077 + test_value(model4.table().elements(), 100, "Test table() for XML constructor");
  2078 + test_value(model4.eval(dir, energy, time), 1.0140199, "Test eval() for XML constructor");
  2079 + test_value(model4.norm(), 1.0, "Test norm() for XML constructor");
  2080 +
  2081 + // Test lookup constructor
  2082 + GEbounds ebds(10, GEnergy(0.2,"TeV"), GEnergy(50.0,"TeV"));
  2083 + GCTAModelSpatialLookup model5(2.0, 0.2, ebds);
  2084 + test_assert(model5.table().has_table("BKG"), "Test table() for lookup constructor");
  2085 + test_value(model5.table().elements(), 100, "Test table() for lookup constructor");
  2086 + test_value(model5.eval(dir, energy, time), 0.0, "Test eval() for lookup constructor");
  2087 + test_value(model5.norm(), 1.0, "Test norm() for lookup constructor");
  2088 +
  2089 + // Test copy constructor
  2090 + GCTAModelSpatialLookup model6(model4);
  2091 + test_assert(model6.table().has_table("BKG"), "Test table() for copy constructor");
  2092 + test_value(model6.table().elements(), 100, "Test table() for copy constructor");
  2093 + test_value(model6.eval(dir, energy, time), 1.0140199, "Test eval() for copy constructor");
  2094 + test_value(model6.norm(), 1.0, "Test norm() for copy constructor");
  2095 +
  2096 + // Test norm() setter and getter
  2097 + model6.norm(0.5);
  2098 + test_value(model6.norm(), 0.5, "Test norm() setter and getter");
  2099 +
  2100 + // Return
  2101 + return;
  2102 +}
  2103 +
  2104 +
  2105 +/***********************************************************************//**
2035 2106 * @brief Test GCTAModelSpatialGaussSpectrum class
2036 2107 ***************************************************************************/
2037 2108 void TestGCTAModel::test_spatial_gauss_spectrum(void)
... ...
inst/cta/test/test_CTA.hpp
  View file @ce238e9
... ... @@ -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_lookup(void);
128 129 void test_spatial_gauss_spectrum(void);
129 130 };
130 131  
... ...