From ce238e9657dc2d5e59af5a3af284557b11268759 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 26 Apr 2019 21:56:36 +0200 Subject: [PATCH] Add GCTAModelSpatialLookup class to handle CTA background lookups (#2878) --- ChangeLog | 3 ++- NEWS | 6 +++++- doc/source/admin/release_history/1.6.rst | 2 ++ inst/cta/Makefile.am | 2 ++ inst/cta/include/GCTALib.hpp | 1 + inst/cta/include/GCTAModelSpatialLookup.hpp | 203 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ inst/cta/pyext/GCTAModelSpatialLookup.i | 101 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ inst/cta/pyext/cta.i | 1 + inst/cta/src/GCTAModelSpatialLookup.cpp | 1069 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ inst/cta/test/data/bkg_lookup.fits | Bin 0 -> 8640 bytes inst/cta/test/data/cta_model_bkg_lookup.xml | 13 +++++++++++++ inst/cta/test/test_CTA.cpp | 73 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++- inst/cta/test/test_CTA.hpp | 1 + 13 files changed, 1472 insertions(+), 3 deletions(-) create mode 100644 inst/cta/include/GCTAModelSpatialLookup.hpp create mode 100644 inst/cta/pyext/GCTAModelSpatialLookup.i create mode 100644 inst/cta/src/GCTAModelSpatialLookup.cpp create mode 100644 inst/cta/test/data/bkg_lookup.fits create mode 100644 inst/cta/test/data/cta_model_bkg_lookup.xml diff --git a/ChangeLog b/ChangeLog index 15d2105..3d25da0 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2019-04-13 +2019-04-26 * Version 1.6.0 released ======================== + Add GCTAModelSpatialLookup class to handle CTA background lookups (#2878) Implement energy-dependent Gaussian CTA spatial background model (#2868) Take instrument scale factors in On/Off analysis into account (#2867) GWcs::wcs_set() method now sets radesys and equinox members (#2866) diff --git a/NEWS b/NEWS index f427bc8..f9c300e 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.6.0 -13 April 2019 +26 April 2019 1. Introduction @@ -21,6 +21,7 @@ The following classes have been added: - GCTAModelSpatialMultiplicative - GCTAModelSpatialRegistry - GCTAModelSpatialGaussSpectrum +- GCTAModelSpatialLookup The following classes have been removed: - GCTACubeSourceDiffuse @@ -332,6 +333,9 @@ unit tests for that response functions were added. 20. CTA interface ----------------- +Add the GCTAModelSpatialLookup class that implements a lookup table for the +energy-dependent spatial modelling of the background (#2878). + Add the GCTAModelSpatialGaussSpectrum class that implements an energy-dependent Gaussian spatial background model (#2868). diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 3d0cd75..53dbec8 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -55,6 +55,8 @@ Bug fixes Improvements ------------ +* [`2878 `_] - + Add ``GCTAModelSpatialLookup`` class to handle CTA background lookups * [`2868 `_] - Implement energy-dependent Gaussian CTA spatial background model * [`2867 `_] - diff --git a/inst/cta/Makefile.am b/inst/cta/Makefile.am index d8f400d..a17abee 100644 --- a/inst/cta/Makefile.am +++ b/inst/cta/Makefile.am @@ -75,6 +75,7 @@ sources = src/GCTASupport.cpp \ src/GCTAPointing.cpp \ src/GCTAModelBackground.cpp \ src/GCTAModelSpatial.cpp \ + src/GCTAModelSpatialLookup.cpp \ src/GCTAModelSpatialGaussSpectrum.cpp \ src/GCTAModelSpatialGradient.cpp \ src/GCTAModelSpatialMultiplicative.cpp \ @@ -130,6 +131,7 @@ pkginclude_HEADERS = include/GCTATypemaps.hpp \ include/GCTACubeSource.hpp \ include/GCTAModelBackground.hpp \ include/GCTAModelSpatial.hpp \ + include/GCTAModelSpatialLookup.hpp \ include/GCTAModelSpatialGaussSpectrum.hpp \ include/GCTAModelSpatialGradient.hpp \ include/GCTAModelSpatialMultiplicative.hpp \ diff --git a/inst/cta/include/GCTALib.hpp b/inst/cta/include/GCTALib.hpp index bef4820..50daafa 100644 --- a/inst/cta/include/GCTALib.hpp +++ b/inst/cta/include/GCTALib.hpp @@ -73,6 +73,7 @@ #include "GCTACubeSource.hpp" #include "GCTAModelBackground.hpp" #include "GCTAModelSpatial.hpp" +#include "GCTAModelSpatialLookup.hpp" #include "GCTAModelSpatialGaussSpectrum.hpp" #include "GCTAModelSpatialGradient.hpp" #include "GCTAModelSpatialMultiplicative.hpp" diff --git a/inst/cta/include/GCTAModelSpatialLookup.hpp b/inst/cta/include/GCTAModelSpatialLookup.hpp new file mode 100644 index 0000000..3462f01 --- /dev/null +++ b/inst/cta/include/GCTAModelSpatialLookup.hpp @@ -0,0 +1,203 @@ +/*************************************************************************** + * GCTAModelSpatialLookup.hpp - Spatial lookup table model * + * ----------------------------------------------------------------------- * + * copyright (C) 2019 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAModelSpatialLookup.hpp + * @brief Spatial lookup table model interface definition + * @author Juergen Knoedlseder + */ + +#ifndef GCTAMODELSPATIALLOOKUP_HPP +#define GCTAMODELSPATIALLOOKUP_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GCTAModelSpatial.hpp" +#include "GCTAResponseTable.hpp" + +/* __ Forward declarations _______________________________________________ */ +class GFilename; +class GXmlElement; +class GCTAInstDir; +class GEnergy; +class GTime; +class GEbounds; +class GFitsTable; +class GFitsBinTable; +class GObservations; +class GCTAObservation; + +/* __ Constants __________________________________________________________ */ +namespace gammalib { + const std::string extname_cta_spatial_lookup = "RADIAL BACKGROUND LOOKUP"; +} + + +/***********************************************************************//** + * @class GCTAModelSpatialLookup + * + * @brief Spatial lookup table model class + ***************************************************************************/ +class GCTAModelSpatialLookup : public GCTAModelSpatial { + +public: + // Constructors and destructors + GCTAModelSpatialLookup(void); + explicit GCTAModelSpatialLookup(const GFilename& filename); + explicit GCTAModelSpatialLookup(const GCTAResponseTable& table); + explicit GCTAModelSpatialLookup(const GXmlElement& xml); + GCTAModelSpatialLookup(const double& maxrad, + const double& radbin, + const GEbounds& ebds); + GCTAModelSpatialLookup(const GCTAModelSpatialLookup& model); + virtual ~GCTAModelSpatialLookup(void); + + // Operators + virtual GCTAModelSpatialLookup& operator=(const GCTAModelSpatialLookup& model); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GCTAModelSpatialLookup* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double eval(const GCTAInstDir& dir, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual double mc_max_value(const GCTAObservation& obs) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + virtual std::string print(const GChatter& chatter = NORMAL) const; + + // Other methods + void fill(const GCTAObservation& obs); + void fill(const GObservations& obs); + const GCTAResponseTable& table(void) const; + void table(const GCTAResponseTable& table); + void read(const GFitsTable& table); + void write(GFitsBinTable& table) const; + void load(const GFilename& filename); + void save(const GFilename& filename, + const bool& clobber = false) const; + double norm(void) const; + void norm(const double& norm); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GCTAModelSpatialLookup& model); + void free_members(void); + void prepare_table(void); + void normalise_table(void); + int table_index(const int& ienergy, const int& itheta) const; + void fill_buffer(const GCTAObservation& obs, std::vector& buffer); + void set_from_buffer(const std::vector& buffer); + + // Protected members + GFilename m_filename; //!< Name of lookup table + GCTAResponseTable m_lookup; //!< Lookup table + GModelPar m_norm; //!< Normalization factor + int m_inx_energy; //!< Energy index + int m_inx_theta; //!< Theta index +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GCTAModelSpatialLookup"). + ***************************************************************************/ +inline +std::string GCTAModelSpatialLookup::classname(void) const +{ + return ("GCTAModelSpatialLookup"); +} + + +/***********************************************************************//** + * @brief Return model type + * + * @return Model type "LookupTable". + ***************************************************************************/ +inline +std::string GCTAModelSpatialLookup::type(void) const +{ + return ("LookupTable"); +} + + +/***********************************************************************//** + * @brief Return maximum function value for Monte Carlo simulations + * + * @param[in] obs CTA Observation. + * @return Maximum function value for Monte Carlo simulations. + * + * This method always returns 1. + ***************************************************************************/ +inline +double GCTAModelSpatialLookup::mc_max_value(const GCTAObservation& obs) const +{ + return 1.0; +} + + +/***********************************************************************//** + * @brief Return lookup table + * + * @return Lookup table. + * + * Returns the lookup table. + ***************************************************************************/ +inline +const GCTAResponseTable& GCTAModelSpatialLookup::table(void) const +{ + return m_lookup; +} + + +/***********************************************************************//** + * @brief Get lookup table model normalisation + * + * @return Lookup table model normalisation. + * + * Returns the normalisation of the lookup table model. + ***************************************************************************/ +inline +double GCTAModelSpatialLookup::norm(void) const +{ + return (m_norm.value()); +} + + +/***********************************************************************//** + * @brief Set lookup table model normalisation + * + * @param[in] norm Lookup table model normalisation. + * + * Set the normalisation of the lookup table model. + ***************************************************************************/ +inline +void GCTAModelSpatialLookup::norm(const double& norm) +{ + m_norm.value(norm); + return; +} + +#endif /* GCTAMODELSPATIALLOOKUP_HPP */ diff --git a/inst/cta/pyext/GCTAModelSpatialLookup.i b/inst/cta/pyext/GCTAModelSpatialLookup.i new file mode 100644 index 0000000..cfde38d --- /dev/null +++ b/inst/cta/pyext/GCTAModelSpatialLookup.i @@ -0,0 +1,101 @@ +/*************************************************************************** + * GCTAModelSpatialLookup.i - Spatial lookup table model * + * ----------------------------------------------------------------------- * + * copyright (C) 2019 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAModelSpatialLookup.i + * @brief Spatial lookup table model interface definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GCTAModelSpatialLookup.hpp" +%} + + +/***********************************************************************//** + * @class GCTAModelSpatialLookup + * + * @brief Spatial lookup table model class + ***************************************************************************/ +class GCTAModelSpatialLookup : public GCTAModelSpatial { + +public: + // Constructors and destructors + GCTAModelSpatialLookup(void); + explicit GCTAModelSpatialLookup(const GFilename& filename); + explicit GCTAModelSpatialLookup(const GCTAResponseTable& table); + explicit GCTAModelSpatialLookup(const GXmlElement& xml); + GCTAModelSpatialLookup(const double& radmax, + const double& radbin, + const GEbounds& ebds); + GCTAModelSpatialLookup(const GCTAModelSpatialLookup& model); + virtual ~GCTAModelSpatialLookup(void); + + // Operators + virtual GCTAModelSpatialLookup& operator=(const GCTAModelSpatialLookup& model); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GCTAModelSpatialLookup* clone(void) const; + virtual std::string classname(void) const; + virtual std::string type(void) const; + virtual double eval(const GCTAInstDir& dir, + const GEnergy& energy, + const GTime& time, + const bool& gradients = false) const; + virtual double mc_max_value(const GCTAObservation& obs) const; + virtual void read(const GXmlElement& xml); + virtual void write(GXmlElement& xml) const; + + // Other methods + void fill(const GCTAObservation& obs); + void fill(const GObservations& obs); + const GCTAResponseTable& table(void) const; + void table(const GCTAResponseTable& table); + void read(const GFitsTable& table); + void write(GFitsBinTable& table) const; + void load(const GFilename& filename); + void save(const GFilename& filename, + const bool& clobber = false) const; + double norm(void) const; + void norm(const double& norm); +}; + + +/***********************************************************************//** + * @brief GCTAModelSpatialLookup class extension + ***************************************************************************/ +%extend GCTAModelSpatialLookup { + GCTAModelSpatialLookup copy() { + return (*self); + } +%pythoncode { + def __getstate__(self): + xml = gammalib.GXmlElement() + self.write(xml) + state = (xml,) + return state + def __setstate__(self, state): + if state[0].elements('parameter') == 0: + self.__init__() + else: + self.__init__(state[0]) +} +}; diff --git a/inst/cta/pyext/cta.i b/inst/cta/pyext/cta.i index d2d4a7b..750bb5f 100644 --- a/inst/cta/pyext/cta.i +++ b/inst/cta/pyext/cta.i @@ -308,6 +308,7 @@ %include "GCTAModelIrfBackground.i" %include "GCTAModelAeffBackground.i" %include "GCTAModelSpatial.i" +%include "GCTAModelSpatialLookup.i" %include "GCTAModelSpatialGaussSpectrum.i" %include "GCTAModelSpatialGradient.i" %include "GCTAModelSpatialMultiplicative.i" diff --git a/inst/cta/src/GCTAModelSpatialLookup.cpp b/inst/cta/src/GCTAModelSpatialLookup.cpp new file mode 100644 index 0000000..ca08148 --- /dev/null +++ b/inst/cta/src/GCTAModelSpatialLookup.cpp @@ -0,0 +1,1069 @@ +/*************************************************************************** + * GCTAModelSpatialLookup.cpp - Spatial lookup table model * + * ----------------------------------------------------------------------- * + * copyright (C) 2019 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file GCTAModelSpatialLookup.cpp + * @brief Spatial lookup table model interface implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "GTools.hpp" +#include "GMath.hpp" +#include "GEnergy.hpp" +#include "GTime.hpp" +#include "GEbounds.hpp" +#include "GFits.hpp" +#include "GFitsTable.hpp" +#include "GFitsBinTable.hpp" +#include "GObservations.hpp" +#include "GCTAObservation.hpp" +#include "GCTAInstDir.hpp" +#include "GCTAModelSpatialLookup.hpp" +#include "GCTAModelSpatialRegistry.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GCTAModelSpatialLookup g_cta_spatial_lookup_seed; +const GCTAModelSpatialRegistry g_cta_spatial_lookup_registry(&g_cta_spatial_lookup_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_CONSTRUCTOR "GCTAModelSpatialLookup(double&, double&, GEbounds&)" +#define G_READ_XML "GCTAModelSpatialLookup::read(GXmlElement&)" +#define G_WRITE_XML "GCTAModelSpatialLookup::write(GXmlElement&)" +#define G_FILL_1 "GCTAModelSpatialLookup::fill(GCTAObservation&)" +#define G_FILL_2 "GCTAModelSpatialLookup::fill(GObservations&)" +#define G_READ_TABLE "GCTAModelSpatialLookup::read(GFitsTable&)" +#define G_LOAD "GCTAModelSpatialLookup::load(GFilename&)" +#define G_FILL_BUFFER "GCTAModelSpatialLookup::fill_buffer(GCTAObservation&,"\ + " std::vector&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GCTAModelSpatialLookup::GCTAModelSpatialLookup(void) : GCTAModelSpatial() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Lookup table file constructor + * + * @param[in] filename Lookup table file name. + * + * Creates instance of a spatial lookup table model from a lookup table FITS + * file. + ***************************************************************************/ +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const GFilename& filename) : + GCTAModelSpatial() +{ + // Initialise members + init_members(); + + // Load lookup table + load(filename); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Response table constructor + * + * @param[in] table Response table. + * + * Creates instance of a spatial lookup table model from a response table. + ***************************************************************************/ +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const GCTAResponseTable& table) : + GCTAModelSpatial() +{ + // Initialise members + init_members(); + + // Set lookup table + this->table(table); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Creates instance of a spatial lookup table model from an XML element. See + * the read() method for information about the expected structure of the XML + * element. + ***************************************************************************/ +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const GXmlElement& xml) : + GCTAModelSpatial() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Lookup table constructor + * + * @param[in] maxrad Maximum radius (deg). + * @param[in] radbin Radial binsize (deg). + * @param[in] ebds Energy boundaries. + * + * @exception GException::invalid_argument + * Non positive radial binsize specified. + * + * Creates empty instance of a spatial lookup table model using the maximum + * lookup radius @p maxrad, the radial lookup binsize @p radbin, and some + * energy boundaries @p ebds. Following construction, the instance may be + * filled using the fill() methods. + * + * The method throws an exception if a non-positive radial bin size is + * specified. + ***************************************************************************/ +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const double& maxrad, + const double& radbin, + const GEbounds& ebds) : + GCTAModelSpatial() +{ + // Throw an exception if the radian binsize is not positive + if (radbin <= 0.0) { + std::string msg = "Non-positive radial bin size "+ + gammalib::str(radbin)+" specified."; + throw GException::invalid_argument(G_CONSTRUCTOR, msg); + } + + // Initialise members + init_members(); + + // Setup energy axis vectors + std::vector eng_lo(ebds.size()); + std::vector eng_hi(ebds.size()); + for (int i = 0; i < ebds.size(); ++i) { + eng_lo[i] = ebds.emin(i).TeV(); + eng_hi[i] = ebds.emax(i).TeV(); + } + + // Setup radial axis vector + int nrad = maxrad/radbin; + std::vector rad_lo(nrad); + std::vector rad_hi(nrad); + for (int i = 0; i < nrad; ++i) { + rad_lo[i] = double(i)*radbin; + rad_hi[i] = double(i+1)*radbin; + } + + // Setup response table + m_lookup.append_axis(eng_lo, eng_hi, "ENERG", "TeV"); + m_lookup.append_axis(rad_lo, rad_hi, "THETA", "deg"); + m_lookup.append_table("BKG", "1/sr"); + + // Prepare table + prepare_table(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Lookup table spatial model. + ***************************************************************************/ +GCTAModelSpatialLookup::GCTAModelSpatialLookup(const GCTAModelSpatialLookup& model) : + GCTAModelSpatial(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GCTAModelSpatialLookup::~GCTAModelSpatialLookup(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Lookup table spatial model. + ***************************************************************************/ +GCTAModelSpatialLookup& GCTAModelSpatialLookup::operator=(const GCTAModelSpatialLookup& model) +{ + // Execute only if object is not identical + if (this != &model) { + + // Copy base class members + this->GCTAModelSpatial::operator=(model); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + } // endif: object was not identical + + // Return + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear instance + ***************************************************************************/ +void GCTAModelSpatialLookup::clear(void) +{ + // Free class members (base and derived classes, derived class first) + free_members(); + this->GCTAModelSpatial::free_members(); + + // Initialise members + this->GCTAModelSpatial::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Clone instance + ***************************************************************************/ +GCTAModelSpatialLookup* GCTAModelSpatialLookup::clone(void) const +{ + return new GCTAModelSpatialLookup(*this); +} + + +/***********************************************************************//** + * @brief Evaluate function + * + * @param[in] dir Event direction. + * @param[in] energy Event energy. + * @param[in] time Event time. + * @param[in] gradients Compute gradients? + * @return Function value + * + * Evaluates the lookup table model for a given event direction and energy. + * The evaluation is done using bi-linear interpolation in the logarithm + * of the energy and the offset angle radius. The lookup table value is + * multiplied with a normalisation parameter. The method always returns a + * non-negative result. + * + * If the @p gradients flag is true the method will also compute the partial + * derivatives of the normalisation parameter. + ***************************************************************************/ +double GCTAModelSpatialLookup::eval(const GCTAInstDir& dir, + const GEnergy& energy, + const GTime& time, + const bool& gradients) const +{ + // Compute offset angle in radians + double offset = dir.theta(); + + // Determine function value by bilinear interpolation + double value = ((m_lookup.elements() > 0) && + (m_lookup.tables() > 0)) + ? m_lookup(0, energy.log10TeV(), offset) : 0.0; + + // Make sure that function does not become negative + if (value < 0.0) { + value = 0.0; + } + + // Optionally compute partial derivatives + if (gradients) { + + // Compute partial derivative of the normalisation parameter + double g_norm = (m_norm.is_free()) ? value * m_norm.scale() : 0.0; + + // Set gradient + m_norm.factor_gradient(g_norm); + + } // endif: computed partial derivatives + + // Return intensity times normalization factor + return (value * m_norm.value()); +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * Read the lookup table model information from an XML element. The XML + * element needs to be of the following format: + * + * + * + * + * + * The @p file attribute provides the filename of the lookup table FITS file. + * The filename may be either an absolute filename (starting with '/') or a + * relative filename. If no access path is given, the file is expected to + * reside in the same location as the XML file. + ***************************************************************************/ +void GCTAModelSpatialLookup::read(const GXmlElement& xml) +{ + // Clear model + clear(); + + // Load lookup table + load(gammalib::xml_file_expand(xml, xml.attribute("file"))); + + // Get parameter pointers + const GXmlElement* norm = gammalib::xml_get_par(G_READ_XML, xml, m_norm.name()); + + // Read parameters + m_norm.read(*norm); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write model into XML element + * + * @param[in] xml XML element. + * + * @exception GException::invalid_value + * Spatial model is not of valid type. + * + * Write the lookup table model information into an XML element. The XML + * element will be of the following format: + * + * + * + * + * + * The @p file attribute provides the filename of the lookup table FITS file. + * The filename may be either an absolute filename (starting with '/') or a + * relative filename. If no access path is given, the file is expected to + * reside in the same location as the XML file. + ***************************************************************************/ +void GCTAModelSpatialLookup::write(GXmlElement& xml) const +{ + // Set model type + if (xml.attribute("type") == "") { + xml.attribute("type", type()); + } + + // Verify model type + if (xml.attribute("type") != type()) { + std::string msg = "Spatial model \""+xml.attribute("type")+ + "\" is not of type \""+type()+"\"."; + throw GException::invalid_value(G_WRITE_XML, msg); + } + + // Set lookup table file name + xml.attribute("file", gammalib::xml_file_reduce(xml, m_filename)); + + // Get XML parameters + GXmlElement* norm = gammalib::xml_need_par(G_WRITE_XML, xml, m_norm.name()); + + // Write parameters + m_norm.write(*norm); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Fill lookup table with events from one CTA observation + * + * @param[in] obs CTA observation. + * + * @exception GException::invalid_value + * No lookup table has been defined. + * + * Fills lookup table with events from one CTA observation. The method needs + * that a lookup table is defined before. In case that no lookup table is + * defined the method throws an exception. + ***************************************************************************/ +void GCTAModelSpatialLookup::fill(const GCTAObservation& obs) +{ + // Throw an exception if response table is not initialised + if (m_lookup.elements() == 0) { + std::string msg = "Attempting to fill a yet undefined lookup table. " + "Please allocate a valid lookup table before " + "calling this method."; + throw GException::invalid_value(G_FILL_1, msg); + } + if (m_lookup.tables() == 0) { + std::string msg = "Attempting to fill non-existing lookup table. " + "While the lookup table axes were defined, no table " + "exists. Please allocate a valid lookup table before " + "calling this method."; + throw GException::invalid_value(G_FILL_1, msg); + } + + // Allocate buffer for events + std::vector counts(m_lookup.elements()); + + // Fill counts into buffer + fill_buffer(obs, counts); + + // Set lookup table from buffer + set_from_buffer(counts); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Fill lookup table with CTA events in observation container + * + * @param[in] obs Observation container. + * + * @exception GException::invalid_value + * No lookup table has been defined. + * + * Fills lookup table with CTA events in observation container. The method + * needs that a lookup table is defined before. In case that no lookup table + * is defined the method throws an exception. + ***************************************************************************/ +void GCTAModelSpatialLookup::fill(const GObservations& obs) +{ + // Throw an exception if response table is not initialised + if (m_lookup.elements() == 0) { + std::string msg = "Attempting to fill a yet undefined lookup table. " + "Please allocate a valid lookup table before " + "calling this method."; + throw GException::invalid_value(G_FILL_2, msg); + } + if (m_lookup.tables() == 0) { + std::string msg = "Attempting to fill non-existing lookup table. " + "While the lookup table axes were defined, no table " + "exists. Please allocate a valid lookup table before " + "calling this method."; + throw GException::invalid_value(G_FILL_2, msg); + } + + // Allocate buffer for events + std::vector counts(m_lookup.elements()); + + // Loop over all observations in container + for (int i = 0; i < obs.size(); ++i) { + + // Get observation and continue only if it is a CTA observation + const GCTAObservation* cta = dynamic_cast(obs[i]); + + // Skip observation if it's not CTA + if (cta == NULL) { + continue; + } + + // Skip observation if we have a binned observation + if (cta->eventtype() == "CountsCube") { + continue; + } + + // Fill counts into buffer + fill_buffer(*cta, counts); + + } + + // Set lookup table from buffer + set_from_buffer(counts); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Assign lookup table + * + * @param[in] table Lookup table. + * + * Assigns lookup table. + ***************************************************************************/ +void GCTAModelSpatialLookup::table(const GCTAResponseTable& table) +{ + // Assign lookup table + m_lookup = table; + + // Prepare table + prepare_table(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Read lookup table from FITS table + * + * @param[in] table FITS table. + * + * @exception GException::invalid_value + * Response table is not two-dimensional. + * + * Reads the lookup table form the FITS @p table. The following column names + * are mandatory: + * + * ENERG_LO - Energy lower bin boundaries + * ENERG_HI - Energy upper bin boundaries + * THETA_LO - Offset angle lower bin boundaries + * THETA_HI - Offset angle upper bin boundaries + * + * After reading, the method assures that the lookup table is properly + * normalised, i.e. that for each energy vector the radial profile has a + * maximum value of one. + ***************************************************************************/ +void GCTAModelSpatialLookup::read(const GFitsTable& table) +{ + // Clear lookup table + m_lookup.clear(); + + // Read lookup table + m_lookup.read(table); + + // Throw an exception if the table is not two-dimensional + if (m_lookup.axes() != 2) { + std::string msg = "Expected two-dimensional lookup table but found "+ + gammalib::str(m_lookup.axes())+" dimensions. Please " + "specify a two-dimensional lookup table."; + throw GException::invalid_value(G_READ_TABLE, msg); + } + + // Prepare table + prepare_table(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Write lookup table into FITS binary table + * + * @param[in] table FITS binary table. + * + * Writes lookup table into a FITS binary @p table. + ***************************************************************************/ +void GCTAModelSpatialLookup::write(GFitsBinTable& table) const +{ + // Write lookup table + m_lookup.write(table); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Load lookup table + * + * @param[in] filename Load lookup table file name. + * + * Loads lookup table from FITS file. + ***************************************************************************/ +void GCTAModelSpatialLookup::load(const GFilename& filename) +{ + // Clear sigma spectrum + clear(); + + // Continue only if filename is not empty + if (!filename.is_empty()) { + + // Open FITS file + GFits fits(filename); + + // Open FITS table + const GFitsTable* table = NULL; + if (filename.has_extname()) { + table = fits.table(filename.extname()); + } + else if (filename.has_extno()) { + table = fits.table(filename.extno()); + } + else { + table = fits.table(1); + } + + // If no FITS table is found then throw an exception + if (table == NULL) { + std::string msg = "No lookup table found in file \""+ + filename.url()+"\"."; + throw GException::file_error(G_LOAD, msg); + } + + // Read lookup table + read(*table); + + // Close FITS file + fits.close(); + + // Store filename + m_filename = filename; + + } // endif: filename was not empty + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save lookup table into FITS file + * + * @param[in] filename FITS file name. + * @param[in] clobber Overwrite existing file? + * + * Saves lookup table into a FITS file. If a file with the given @p filename + * does not yet exist it will be created, otherwise the method opens the + * existing file. The method will create a (or replace an existing) + * effective area extension. The extension name can be specified as part + * of the @p filename, or if no extension name is given, is assumed to be + * `RADIAL BACKGROUND LOOKUP`. + * + * An existing file will only be modified if the @p clobber flag is set to + * true. + ***************************************************************************/ +void GCTAModelSpatialLookup::save(const GFilename& filename, const bool& clobber) const +{ + // Get extension name + std::string extname = filename.extname(gammalib::extname_cta_spatial_lookup); + + // Open or create FITS file (without extension name since the requested + // extension may not yet exist in the file) + GFits fits(filename.url(), true); + + // Remove extension if it exists already + if (fits.contains(extname)) { + fits.remove(extname); + } + + // Create binary table + GFitsBinTable table; + + // Write the lookup table + write(table); + + // Set binary table extension name + table.extname(extname); + + // Append table to FITS file + fits.append(table); + + // Save to file + fits.save(clobber); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print lookup table information + * + * @param[in] chatter Chattiness. + * @return String containing lookup table information. + ***************************************************************************/ +std::string GCTAModelSpatialLookup::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Initialise information + int nebins = 0; + int nthetabins = 0; + double emin = 0.0; + double emax = 0.0; + double omin = 0.0; + double omax = 0.0; + + // Extract information if there are axes in the response table + if (m_lookup.axes() > 0) { + nebins = m_lookup.axis_bins(m_inx_energy); + nthetabins = m_lookup.axis_bins(m_inx_theta); + emin = m_lookup.axis_lo(m_inx_energy,0); + emax = m_lookup.axis_hi(m_inx_energy,nebins-1); + omin = m_lookup.axis_lo(m_inx_theta,0); + omax = m_lookup.axis_hi(m_inx_theta,nthetabins-1); + } + + // Append header + result.append("=== GCTAModelSpatialLookup ==="); + + // Append information + result.append("\n"+gammalib::parformat("Filename")+m_filename); + result.append("\n"+gammalib::parformat("Number of energy bins") + + gammalib::str(nebins)); + result.append("\n"+gammalib::parformat("Number of offset bins") + + gammalib::str(nthetabins)); + result.append("\n"+gammalib::parformat("Energy range")); + result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV"); + result.append("\n"+gammalib::parformat("Offset angle range")); + result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg"); + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GCTAModelSpatialLookup::init_members(void) +{ + // Initialise members + m_filename.clear(); + m_lookup.clear(); + m_inx_energy = 0; + m_inx_theta = 1; + + // Initialise lookup table normalisation + m_norm.clear(); + m_norm.name("Normalization"); + m_norm.value(1.0); + m_norm.scale(1.0); + m_norm.range(0.001, 1000.0); + m_norm.gradient(0.0); + m_norm.fix(); + m_norm.has_grad(true); + + // Set parameter pointer(s) + m_pars.clear(); + m_pars.push_back(&m_norm); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Lookup table model. + ***************************************************************************/ +void GCTAModelSpatialLookup::copy_members(const GCTAModelSpatialLookup& model) +{ + // Copy members + m_filename = model.m_filename; + m_lookup = model.m_lookup; + m_norm = model.m_norm; + m_inx_energy = model.m_inx_energy; + m_inx_theta = model.m_inx_theta; + + // Set parameter pointer(s) + m_pars.clear(); + m_pars.push_back(&m_norm); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GCTAModelSpatialLookup::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Prepare lookup table indices + * + * Prepares a lookup table for usage. + * + * The method sets the data members m_inx_energy and m_inx_theta which + * determine the axes order. + * + * The method furthermore sets the energy axis to logarithmic scale and the + * offset angles to radians. + * + * Finally, the method assures that the maximum values for each radial + * profile are one. + ***************************************************************************/ +void GCTAModelSpatialLookup::prepare_table(void) +{ + // Get mandatory indices (throw exception if not found) + m_inx_energy = m_lookup.axis("ENERG"); + m_inx_theta = m_lookup.axis("THETA"); + + // Set energy axis to logarithmic scale + m_lookup.axis_log10(m_inx_energy); + + // Set offset angle axis to radians + m_lookup.axis_radians(m_inx_theta); + + // Normalise radial profiles + normalise_table(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Normalise lookup table + * + * Normalises lookup table so that each radial profile has a maximum value + * of one. + ***************************************************************************/ +void GCTAModelSpatialLookup::normalise_table(void) +{ + // Get axes dimensions + int energy_size = m_lookup.axis_bins(m_inx_energy); + int theta_size = m_lookup.axis_bins(m_inx_theta); + + // Loop over energies + for (int i_energy = 0; i_energy < energy_size; ++i_energy) { + + // Initialise maximum radial profile value + double rad_max = 0.0; + + // Determine maximum radial profile value + for (int i_theta = 0; i_theta < theta_size; ++i_theta) { + int i = table_index(i_energy, i_theta); + if (m_lookup(0,i) > rad_max) { + rad_max = m_lookup(0,i); + } + } + + // If maximum value is positive then normalise profile + if (rad_max > 0.0) { + for (int i_theta = 0; i_theta < theta_size; ++i_theta) { + int i = table_index(i_energy, i_theta); + m_lookup(0,i) /= rad_max; + } + } + + } // endfor: looped over all energies + + // Return + return; +} + + +/***********************************************************************//** + * @brief Return index of lookup table element + * + * @param[in] ienergy Energy index. + * @param[in] itheta Offset index. + * @return Index of lookup table element. + ***************************************************************************/ +int GCTAModelSpatialLookup::table_index(const int& ienergy, + const int& itheta) const +{ + // Set index vector + int inx[2]; + inx[m_inx_energy] = ienergy; + inx[m_inx_theta] = itheta; + + // Compute index + int index = inx[0] + inx[1] * m_lookup.axis_bins(0); + + // Return index + return index; +} + + +/***********************************************************************//** + * @brief Fill buffer from events in CTA observation + * + * @param[in] obs CTA observation. + * @param[in,out] buffer Counts buffer. + ***************************************************************************/ +void GCTAModelSpatialLookup::fill_buffer(const GCTAObservation& obs, + std::vector& buffer) +{ + // Make sure that the observation holds a CTA event list. If this is + // not the case then throw an exception. + const GCTAEventList* events = dynamic_cast(obs.events()); + if (events == NULL) { + std::string msg = "CTA Observation does not contain an event " + "list. An event list is needed to fill the " + "lookup table."; + throw GException::invalid_value(G_FILL_BUFFER, msg); + } + + // Get axes dimensions + int energy_size = m_lookup.axis_bins(m_inx_energy); + int theta_size = m_lookup.axis_bins(m_inx_theta); + + // Fill buffer + for (int i = 0; i < events->size(); ++i) { + + // Get event + const GCTAEventAtom* event = (*events)[i]; + + // Determine offset angle in degrees + GCTAInstDir* inst = (GCTAInstDir*)&(event->dir()); + double theta = inst->theta() * gammalib::rad2deg; + + // Determine offset angle bin + int i_theta = 0; + for (; i_theta < theta_size; ++i_theta) { + if ((theta >= m_lookup.axis_lo(m_inx_theta,i_theta)) && + (theta <= m_lookup.axis_hi(m_inx_theta,i_theta))) { + break; + } + } + + // Skip event if no valid offset angle bin was found + if (i_theta >= theta_size) { + continue; + } + + // Determine event energy in TeV + double energy = event->energy().TeV(); + + // Determine energy bin + int i_energy = 0; + for (; i_energy < energy_size; ++i_energy) { + if ((energy >= m_lookup.axis_lo(m_inx_energy,i_energy)) && + (energy <= m_lookup.axis_hi(m_inx_energy,i_energy))) { + break; + } + } + + // Skip event if no valid offset angle bin was found + if (i_energy >= theta_size) { + continue; + } + + // Fill event into bugger + int inx = table_index(i_energy, i_theta); + buffer[inx] += 1.0; + + } // endfor: looped over all events + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set lookup table from buffer + * + * @param[in] buffer Counts buffer. + ***************************************************************************/ +void GCTAModelSpatialLookup::set_from_buffer(const std::vector& buffer) +{ + // Get axes dimensions + int energy_size = m_lookup.axis_bins(m_inx_energy); + int theta_size = m_lookup.axis_bins(m_inx_theta); + + // Loop over offset angle + for (int i_theta = 0; i_theta < theta_size; ++i_theta) { + + // Compute solid angle for offset angle bin + double theta_min = m_lookup.axis_lo(m_inx_theta,i_theta) * gammalib::deg2rad; + double theta_max = m_lookup.axis_hi(m_inx_theta,i_theta) * gammalib::deg2rad; + double area = gammalib::pi * (theta_max*theta_max - theta_min*theta_min); + + // Fill energy vector into lookup table + for (int i_energy = 0; i_energy < energy_size; ++i_energy) { + int inx = table_index(i_energy, i_theta); + m_lookup(0,inx) = buffer[inx] / area; + } + + } // endfor: looped over offset angles + + // Normalise radial profiles + normalise_table(); + + // Return + return; +} diff --git a/inst/cta/test/data/bkg_lookup.fits b/inst/cta/test/data/bkg_lookup.fits new file mode 100644 index 0000000..faaa4eb Binary files /dev/null and b/inst/cta/test/data/bkg_lookup.fits differ diff --git a/inst/cta/test/data/cta_model_bkg_lookup.xml b/inst/cta/test/data/cta_model_bkg_lookup.xml new file mode 100644 index 0000000..e9d391d --- /dev/null +++ b/inst/cta/test/data/cta_model_bkg_lookup.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index 8170940..575480f 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -63,6 +63,8 @@ const std::string cta_bgd_gauss_e_xml = datadir+"/cta_model_bgd_gauss_e.xml"; const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml"; const std::string cta_irf_bgd_xml = datadir+"/cta_model_irf_bgd.xml"; const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml"; +const std::string cta_bkg_lookup = datadir+"/bkg_lookup.fits"; +const std::string cta_bkg_lookup_xml = datadir+"/cta_model_bkg_lookup.xml"; const std::string cta_caldb_king = cta_caldb+"/data/cta/e/bcf/IFAE20120510_50h_King"; const std::string cta_irf_king = "irf_file.fits"; const std::string cta_psf_table = caldbdir+"/psf_table.fits[PSF_2D_TABLE]"; @@ -70,7 +72,7 @@ const std::string cta_perf_table = caldbdir+"/cta_dummy_irf.dat"; const std::string cta_edisp_rmf = caldbdir+"/dc1/rmf.fits"; const std::string cta_edisp_2D = caldbdir+"/edisp_matrix.fits"; const std::string cta_bgd_3D = caldbdir+"/edisp_matrix.fits"; -const std::string cta_modbck_fit = datadir+"/bg_test.fits"; +//const std::string cta_modbck_fit = datadir+"/bg_test.fits"; /* __ Test files for stacked analysis (based on Prod2::South_0.5h) _______ */ const std::string cta_stacked_xml = datadir+"/stacked_obs.xml"; @@ -177,6 +179,8 @@ void TestGCTAModel::set(void) "Test CTA IRF background model"); append(static_cast(&TestGCTAModel::test_model_aeff_bgd), "Test CTA Aeff background model"); + append(static_cast(&TestGCTAModel::test_spatial_lookup), + "Test GCTAModelSpatialLookup class"); append(static_cast(&TestGCTAModel::test_spatial_gauss_spectrum), "Test GCTAModelSpatialGaussSpectrum class"); @@ -2032,6 +2036,73 @@ void TestGCTAModel::test_model_aeff_bgd(void) /***********************************************************************//** + * @brief Test GCTAModelSpatialLookup class + ***************************************************************************/ +void TestGCTAModel::test_spatial_lookup(void) +{ + // Set instrument direction, energy and time + GCTAInstDir dir(0.0, 0.0); + GEnergy energy(1.0, "TeV"); + GTime time; + + // Test void constructor + GCTAModelSpatialLookup model1; + test_assert(!model1.table().has_table("BKG"), "Test table() for void constructor"); + test_value(model1.table().elements(), 0, "Test table() for void constructor"); + test_value(model1.eval(dir, energy, time), 0.0, "Test eval() for void constructor"); + test_value(model1.norm(), 1.0, "Test norm() for void constructor"); + + // Test filename constructor + GCTAModelSpatialLookup model2(cta_bkg_lookup); + test_assert(model2.table().has_table("BKG"), "Test table() for filename constructor"); + test_value(model2.table().elements(), 100, "Test table() for filename constructor"); + test_value(model2.eval(dir, energy, time), 1.0140199, "Test eval() for filename constructor"); + test_value(model2.norm(), 1.0, "Test norm() for filename constructor"); + + // Test response table constructor + GFits fits(cta_bkg_lookup); + const GFitsTable& hdu = *fits.table(1); + GCTAResponseTable table(hdu); + GCTAModelSpatialLookup model3(table); + test_assert(model3.table().has_table("BKG"), "Test table() for response table constructor"); + test_value(model3.table().elements(), 100, "Test table() for response table constructor"); + test_value(model3.eval(dir, energy, time), 1.0140199, "Test eval() for response table constructor"); + test_value(model3.norm(), 1.0, "Test norm() for response table constructor"); + + // Test XML constructor + GXml xml(cta_bkg_lookup_xml); + GXmlElement* element = xml.element("source_library > source[0] > spatialModel"); + GCTAModelSpatialLookup model4(*element); + test_assert(model4.table().has_table("BKG"), "Test table() for XML constructor"); + test_value(model4.table().elements(), 100, "Test table() for XML constructor"); + test_value(model4.eval(dir, energy, time), 1.0140199, "Test eval() for XML constructor"); + test_value(model4.norm(), 1.0, "Test norm() for XML constructor"); + + // Test lookup constructor + GEbounds ebds(10, GEnergy(0.2,"TeV"), GEnergy(50.0,"TeV")); + GCTAModelSpatialLookup model5(2.0, 0.2, ebds); + test_assert(model5.table().has_table("BKG"), "Test table() for lookup constructor"); + test_value(model5.table().elements(), 100, "Test table() for lookup constructor"); + test_value(model5.eval(dir, energy, time), 0.0, "Test eval() for lookup constructor"); + test_value(model5.norm(), 1.0, "Test norm() for lookup constructor"); + + // Test copy constructor + GCTAModelSpatialLookup model6(model4); + test_assert(model6.table().has_table("BKG"), "Test table() for copy constructor"); + test_value(model6.table().elements(), 100, "Test table() for copy constructor"); + test_value(model6.eval(dir, energy, time), 1.0140199, "Test eval() for copy constructor"); + test_value(model6.norm(), 1.0, "Test norm() for copy constructor"); + + // Test norm() setter and getter + model6.norm(0.5); + test_value(model6.norm(), 0.5, "Test norm() setter and getter"); + + // Return + return; +} + + +/***********************************************************************//** * @brief Test GCTAModelSpatialGaussSpectrum class ***************************************************************************/ void TestGCTAModel::test_spatial_gauss_spectrum(void) diff --git a/inst/cta/test/test_CTA.hpp b/inst/cta/test/test_CTA.hpp index 2e6cd5c..8ae7b3f 100644 --- a/inst/cta/test/test_CTA.hpp +++ b/inst/cta/test/test_CTA.hpp @@ -125,6 +125,7 @@ public: void test_model_cube_bgd(void); void test_model_irf_bgd(void); void test_model_aeff_bgd(void); + void test_spatial_lookup(void); void test_spatial_gauss_spectrum(void); }; -- libgit2 0.21.4