diff --git a/ChangeLog b/ChangeLog index 5c11208..15d2105 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2019-04-09 +2019-04-13 * Version 1.6.0 released ======================== + 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) Add GCTAPointing::is_valid() method and remove table methods diff --git a/NEWS b/NEWS index 1a26ee1..f427bc8 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.6.0 -9 April 2019 +13 April 2019 1. Introduction @@ -20,6 +20,7 @@ The following classes have been added: - GCTAModelSpatial - GCTAModelSpatialMultiplicative - GCTAModelSpatialRegistry +- GCTAModelSpatialGaussSpectrum The following classes have been removed: - GCTACubeSourceDiffuse @@ -331,6 +332,9 @@ unit tests for that response functions were added. 20. CTA interface ----------------- +Add the GCTAModelSpatialGaussSpectrum class that implements an energy-dependent +Gaussian spatial background model (#2868). + Instrument scale factors are now also taken in an On/Off analysis into account. The scale factors were added to the computation in GCTAOnOffObservation::N_gamma() (#2867). diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 7717be4..3d0cd75 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 ------------ +* [`2868 `_] - + Implement energy-dependent Gaussian CTA spatial background model * [`2867 `_] - Take instrument scale factors in On/Off analysis into account * [`2850 `_] - diff --git a/inst/cta/Makefile.am b/inst/cta/Makefile.am index 35696cc..d8f400d 100644 --- a/inst/cta/Makefile.am +++ b/inst/cta/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile.am for cta module # # ------------------------------------------------------------------------- # -# Copyright (C) 2016-2018 Juergen Knoedlseder # +# Copyright (C) 2016-2019 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -75,6 +75,7 @@ sources = src/GCTASupport.cpp \ src/GCTAPointing.cpp \ src/GCTAModelBackground.cpp \ src/GCTAModelSpatial.cpp \ + src/GCTAModelSpatialGaussSpectrum.cpp \ src/GCTAModelSpatialGradient.cpp \ src/GCTAModelSpatialMultiplicative.cpp \ src/GCTAModelSpatialRegistry.cpp \ @@ -129,6 +130,7 @@ pkginclude_HEADERS = include/GCTATypemaps.hpp \ include/GCTACubeSource.hpp \ include/GCTAModelBackground.hpp \ include/GCTAModelSpatial.hpp \ + include/GCTAModelSpatialGaussSpectrum.hpp \ include/GCTAModelSpatialGradient.hpp \ include/GCTAModelSpatialMultiplicative.hpp \ include/GCTAModelSpatialRegistry.hpp \ diff --git a/inst/cta/include/GCTALib.hpp b/inst/cta/include/GCTALib.hpp index 11f4087..bef4820 100644 --- a/inst/cta/include/GCTALib.hpp +++ b/inst/cta/include/GCTALib.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTALib.hpp - CTA Support Header files * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2018 by Juergen Knoedlseder * + * copyright (C) 2010-2019 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -73,6 +73,7 @@ #include "GCTACubeSource.hpp" #include "GCTAModelBackground.hpp" #include "GCTAModelSpatial.hpp" +#include "GCTAModelSpatialGaussSpectrum.hpp" #include "GCTAModelSpatialGradient.hpp" #include "GCTAModelSpatialMultiplicative.hpp" #include "GCTAModelSpatialRegistry.hpp" diff --git a/inst/cta/include/GCTAModelSpatialGaussSpectrum.hpp b/inst/cta/include/GCTAModelSpatialGaussSpectrum.hpp new file mode 100644 index 0000000..a611a1a --- /dev/null +++ b/inst/cta/include/GCTAModelSpatialGaussSpectrum.hpp @@ -0,0 +1,144 @@ +/*************************************************************************** + * GCTAModelSpatialGaussSpectrum.hpp - Spatial energy dependent Gaussian * + * ----------------------------------------------------------------------- * + * 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 GCTAModelSpatialGaussSpectrum.hpp + * @brief Spatial energy dependent Gaussian interface definition + * @author Juergen Knoedlseder + */ + +#ifndef GCTAMODELSPATIALGAUSSSPECTRUM_HPP +#define GCTAMODELSPATIALGAUSSSPECTRUM_HPP + +/* __ Includes ___________________________________________________________ */ +#include +#include "GModelSpectral.hpp" +#include "GCTAModelSpatial.hpp" + +/* __ Forward declarations _______________________________________________ */ +class GXmlElement; +class GCTAInstDir; +class GEnergy; +class GTime; +class GCTAObservation; + + +/***********************************************************************//** + * @class GCTAModelSpatialGaussSpectrum + * + * @brief Spatial energy dependent Gaussian model class + ***************************************************************************/ +class GCTAModelSpatialGaussSpectrum : public GCTAModelSpatial { + +public: + // Constructors and destructors + GCTAModelSpatialGaussSpectrum(void); + explicit GCTAModelSpatialGaussSpectrum(const double& sigma); + explicit GCTAModelSpatialGaussSpectrum(const GModelSpectral& sigma); + explicit GCTAModelSpatialGaussSpectrum(const GXmlElement& xml); + GCTAModelSpatialGaussSpectrum(const GCTAModelSpatialGaussSpectrum& model); + virtual ~GCTAModelSpatialGaussSpectrum(void); + + // Operators + virtual GCTAModelSpatialGaussSpectrum& operator=(const GCTAModelSpatialGaussSpectrum& model); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GCTAModelSpatialGaussSpectrum* 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 + const GModelSpectral* sigma(void) const; + void sigma(const double& sigma); + void sigma(const GModelSpectral& sigma); + +protected: + // Protected methods + void init_members(void); + void copy_members(const GCTAModelSpatialGaussSpectrum& model); + void free_members(void); + void set_pointers(void); + + + // Protected members + GModelSpectral* m_sigma; // Pointer to sigma spectrum +}; + + +/***********************************************************************//** + * @brief Return class name + * + * @return String containing the class name ("GCTAModelSpatialGaussSpectrum"). + ***************************************************************************/ +inline +std::string GCTAModelSpatialGaussSpectrum::classname(void) const +{ + return ("GCTAModelSpatialGaussSpectrum"); +} + + +/***********************************************************************//** + * @brief Return model type + * + * @return Model type "EnergyDependentGaussian". + ***************************************************************************/ +inline +std::string GCTAModelSpatialGaussSpectrum::type(void) const +{ + return ("EnergyDependentGaussian"); +} + + +/***********************************************************************//** + * @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 GCTAModelSpatialGaussSpectrum::mc_max_value(const GCTAObservation& obs) const +{ + return 1.0; +} + + +/***********************************************************************//** + * @brief Return pointer to sigma spectrum + * + * @return Pointer to sigma spectrum. + ***************************************************************************/ +inline +const GModelSpectral* GCTAModelSpatialGaussSpectrum::sigma(void) const +{ + return (m_sigma); +} + +#endif /* GCTAMODELSPATIALGAUSSSPECTRUM_HPP */ diff --git a/inst/cta/pyext/GCTAModelSpatialGaussSpectrum.i b/inst/cta/pyext/GCTAModelSpatialGaussSpectrum.i new file mode 100644 index 0000000..aabe6d1 --- /dev/null +++ b/inst/cta/pyext/GCTAModelSpatialGaussSpectrum.i @@ -0,0 +1,90 @@ +/*************************************************************************** + * GCTAModelSpatialGaussSpectrum.i - Spatial energy dependent Gaussian * + * ----------------------------------------------------------------------- * + * 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 GCTAModelSpatialGaussSpectrum.i + * @brief Spatial energy dependent Gaussian interface definition + * @author Juergen Knoedlseder + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "GCTAModelSpatialGaussSpectrum.hpp" +%} + + +/***********************************************************************//** + * @class GCTAModelSpatialGaussSpectrum + * + * @brief Spatial energy dependent Gaussian model class + ***************************************************************************/ +class GCTAModelSpatialGaussSpectrum : public GCTAModelSpatial { + +public: + // Constructors and destructors + GCTAModelSpatialGaussSpectrum(void); + explicit GCTAModelSpatialGaussSpectrum(const double& sigma); + explicit GCTAModelSpatialGaussSpectrum(const GModelSpectral& sigma); + explicit GCTAModelSpatialGaussSpectrum(const GXmlElement& xml); + GCTAModelSpatialGaussSpectrum(const GCTAModelSpatialGaussSpectrum& model); + virtual ~GCTAModelSpatialGaussSpectrum(void); + + // Operators + virtual GCTAModelSpatialGaussSpectrum& operator=(const GCTAModelSpatialGaussSpectrum& model); + + // Implemented pure virtual methods + virtual void clear(void); + virtual GCTAModelSpatialGaussSpectrum* 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 + const GModelSpectral* sigma(void) const; + void sigma(const double& sigma); + void sigma(const GModelSpectral& sigma); +}; + + +/***********************************************************************//** + * @brief GCTAModelSpatialGaussSpectrum class extension + ***************************************************************************/ +%extend GCTAModelSpatialGaussSpectrum { + GCTAModelSpatialGaussSpectrum 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 2945860..d2d4a7b 100644 --- a/inst/cta/pyext/cta.i +++ b/inst/cta/pyext/cta.i @@ -1,7 +1,7 @@ /*************************************************************************** * cta.i - Cherenkov Telescope Array module * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2018 by Juergen Knoedlseder * + * copyright (C) 2011-2019 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -308,6 +308,7 @@ %include "GCTAModelIrfBackground.i" %include "GCTAModelAeffBackground.i" %include "GCTAModelSpatial.i" +%include "GCTAModelSpatialGaussSpectrum.i" %include "GCTAModelSpatialGradient.i" %include "GCTAModelSpatialMultiplicative.i" %include "GCTAModelSpatialRegistry.i" diff --git a/inst/cta/src/GCTAModelRadialGauss.cpp b/inst/cta/src/GCTAModelRadialGauss.cpp index 187af54..851152b 100644 --- a/inst/cta/src/GCTAModelRadialGauss.cpp +++ b/inst/cta/src/GCTAModelRadialGauss.cpp @@ -533,7 +533,7 @@ void GCTAModelRadialGauss::write(GXmlElement& xml) const /***********************************************************************//** * @brief Print point source information * - * @param[in] chatter Chattiness (defaults to NORMAL). + * @param[in] chatter Chattiness. * @return String containing point source information. ***************************************************************************/ std::string GCTAModelRadialGauss::print(const GChatter& chatter) const diff --git a/inst/cta/src/GCTAModelSpatialGaussSpectrum.cpp b/inst/cta/src/GCTAModelSpatialGaussSpectrum.cpp new file mode 100644 index 0000000..a507c40 --- /dev/null +++ b/inst/cta/src/GCTAModelSpatialGaussSpectrum.cpp @@ -0,0 +1,584 @@ +/*************************************************************************** + * GCTAModelSpatialGaussSpectrum.cpp - Spatial energy dependent Gaussian * + * ----------------------------------------------------------------------- * + * 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 GCTAModelSpatialGaussSpectrum.cpp + * @brief Spatial energy dependent Gaussian interface implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "GTools.hpp" +#include "GMath.hpp" +#include "GCTAInstDir.hpp" +#include "GEnergy.hpp" +#include "GTime.hpp" +#include "GCTAModelSpatialGaussSpectrum.hpp" +#include "GCTAModelSpatialRegistry.hpp" +#include "GModelSpectralRegistry.hpp" +#include "GModelSpectralConst.hpp" + +/* __ Constants __________________________________________________________ */ + +/* __ Globals ____________________________________________________________ */ +const GCTAModelSpatialGaussSpectrum g_cta_spatial_gauss_spec_seed; +const GCTAModelSpatialRegistry g_cta_spatial_gauss_spec_registry(&g_cta_spatial_gauss_spec_seed); + +/* __ Method name definitions ____________________________________________ */ +#define G_READ "GCTAModelSpatialGaussSpectrum::read(GXmlElement&)" +#define G_WRITE "GCTAModelSpatialGaussSpectrum::write(GXmlElement&)" + +/* __ Macros _____________________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + +/* __ Debug definitions __________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(void) : + GCTAModelSpatial() +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Energy-independent Gaussian constructor + * + * @param[in] sigma Gaussian width (degrees\f$^2\f$). + ***************************************************************************/ +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(const double& sigma) : + GCTAModelSpatial() +{ + // Initialise members + init_members(); + + // Assign sigma + this->sigma(sigma); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Gaussian spectrum constructor + * + * @param[in] sigma Spectrum defining the energy-deependent Gaussian width + * (degrees\f$^2\f$). + ***************************************************************************/ +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(const GModelSpectral& sigma) : + GCTAModelSpatial() +{ + // Initialise members + init_members(); + + // Assign sigma spectrum + this->sigma(sigma); + + // Return + return; +} + + +/***********************************************************************//** + * @brief XML constructor + * + * @param[in] xml XML element. + * + * Creates instance of a spatial energy-dependent Gaussian model by + * extracting information from an XML element. See + * GCTAModelSpatialGaussSpectrum::read() for more information about the + * expected structure of the XML element. + ***************************************************************************/ +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(const GXmlElement& xml) : + GCTAModelSpatial() +{ + // Initialise members + init_members(); + + // Read information from XML element + read(xml); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] model Energy-dependent Gaussian model. + ***************************************************************************/ +GCTAModelSpatialGaussSpectrum::GCTAModelSpatialGaussSpectrum(const GCTAModelSpatialGaussSpectrum& model) : + GCTAModelSpatial(model) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(model); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +GCTAModelSpatialGaussSpectrum::~GCTAModelSpatialGaussSpectrum(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] model Energy-dependent Gaussian model. + ***************************************************************************/ +GCTAModelSpatialGaussSpectrum& GCTAModelSpatialGaussSpectrum::operator=(const GCTAModelSpatialGaussSpectrum& 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 GCTAModelSpatialGaussSpectrum::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 + ***************************************************************************/ +GCTAModelSpatialGaussSpectrum* GCTAModelSpatialGaussSpectrum::clone(void) const +{ + return new GCTAModelSpatialGaussSpectrum(*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 energy-dependent Gaussian model for a given event direction. + * The Gaussian model is defined by + * + * \f[f(\theta,E) = \exp \left(-\frac{1}{2} + * \left( \frac{\theta^2}{\sigma(E)} \right)^2 \right)\f] + * + * where + * \f$\theta\f$ is the offset angle (in degrees), and + * \f$\sigma(E)\f$ is the energy-dependent Gaussian width (in degrees\f$^2\f$). + * + * If the @p gradients flag is true the method will also compute the partial + * derivatives of the parameters. The partial derivative of the Gaussian width + * is given by + * + * \f[\frac{df}{d\sigma_v} = f(\theta) \frac{\theta^4}{\sigma^3} \sigma_s\f] + * + * where + * \f$\sigma_v\f$ is the value part, + * \f$\sigma_s\f$ is the scaling part, and + * \f$\sigma = \sigma_v \sigma_s\f$. + * + * Note that this method implements a function which is unity for + * \f$\theta=0\f$. + ***************************************************************************/ +double GCTAModelSpatialGaussSpectrum::eval(const GCTAInstDir& dir, + const GEnergy& energy, + const GTime& time, + const bool& gradients) const +{ + // Initialise value + double value = 0.0; + + // Continue only if a sigma spectrum exists + if (m_sigma != NULL) { + + // Compute offset angle in degrees + double offset = dir.theta() * gammalib::rad2deg; + + // Compute energy dependent sigma value in degrees^2 + double sigma = m_sigma->eval(energy, time, gradients); + + // Compute Gaussian value + double arg = offset * offset / sigma; + double arg2 = arg * arg; + value = std::exp(-0.5 * arg2); + + // Optionally compute partial derivatives + if (gradients) { + + // Compute partial derivatives of the sigma parameter + double g_sigma = value * arg2 / sigma; + + // Loop over model parameters + for (int i = 0; i < m_sigma->size(); ++i) { + + // Get reference to model parameter + GModelPar& par = m_sigma->operator[](i); + + // Scale parameter gradient + par.gradient(par.gradient()*g_sigma); + + } // endfor: loop over model parameters + + } // endif: computed partial derivatives + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GCTAModelSpatialGaussSpectrum::eval"; + std::cout << "(offset=" << offset << "): NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", energy=" << energy; + std::cout << ", sigma=" << sigma; + std::cout << ")" << std::endl; + } + #endif + + } // endif: sigma spectrum existed + + // Return value + return value; +} + + +/***********************************************************************//** + * @brief Read model from XML element + * + * @param[in] xml XML element. + * + * Read the energy-dependent Gaussian spatial model information from an XML + * element. The XML element needs to be of the following format: + * + * + * + * ... + * + * + * + * where any spectral model type can be provided for the @a sigma tag. + ***************************************************************************/ +void GCTAModelSpatialGaussSpectrum::read(const GXmlElement& xml) +{ + // Clear model + clear(); + + // Get pointers on sigma model components + const GXmlElement* sigma = xml.element("sigma", 0); + + // Get sigma model + GModelSpectralRegistry registry; + m_sigma = registry.alloc(*sigma); + + // Set pointers + set_pointers(); + + // 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 energy-dependent Gaussian spatial model information into an XML + * element. The XML element will be of the following format: + * + * + * + * ... + * + * + * + * where the spectral sigma model will be written into the @a sigma tag. + ***************************************************************************/ +void GCTAModelSpatialGaussSpectrum::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, msg); + } + + // Write sigma spectrum if it exists + if (m_sigma != NULL) { + + // Create new sigma spectrum node + xml.append(GXmlElement("sigma")); + + // Get new sigma spectrum node + GXmlElement* sigma = xml.element("sigma", xml.elements("sigma")-1); + + // Write sigma spectrum + m_sigma->write(*sigma); + + } // endif: sigma spectrum was not NULL + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set sigma spectrum from value + * + * @param[in] sigma Constant sigma value. + ***************************************************************************/ +void GCTAModelSpatialGaussSpectrum::sigma(const double& sigma) +{ + // Clear sigma spectrum + clear(); + + // Set constant spectrum + m_sigma = new GModelSpectralConst(sigma); + + // Set parameter pointers + set_pointers(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set sigma spectrum + * + * @param[in] sigma Sigma spectrum. + ***************************************************************************/ +void GCTAModelSpatialGaussSpectrum::sigma(const GModelSpectral& sigma) +{ + // Clear sigma spectrum + clear(); + + // Set spectrum + m_sigma = sigma.clone(); + + // Set parameter pointers + set_pointers(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Print point source information + * + * @param[in] chatter Chattiness. + * @return String containing point source information. + ***************************************************************************/ +std::string GCTAModelSpatialGaussSpectrum::print(const GChatter& chatter) const +{ + // Initialise result string + std::string result; + + // Continue only if chatter is not silent + if (chatter != SILENT) { + + // Append header + result.append("=== GCTAModelSpatialGaussSpectrum ==="); + + // Determine number of sigma parameters + int n_sigma = (m_sigma != NULL) ? m_sigma->size() : 0; + + // Append sigma spectrum type + result.append("\n"+gammalib::parformat("Sigma spectrum")); + if (n_sigma > 0) { + result.append("\""+m_sigma->type()+"\""); + } + + // Append sigma spectrum parameters + for (int i = 0; i < n_sigma; ++i) { + result.append("\n"+(m_sigma)[i].print()); + } + + } // endif: chatter was not silent + + // Return result + return result; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void GCTAModelSpatialGaussSpectrum::init_members(void) +{ + // Initialise sigma spectrum pointer + m_sigma = NULL; + + // Clear parameter pointer(s) + m_pars.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] model Energy-dependent Gaussian model. + ***************************************************************************/ +void GCTAModelSpatialGaussSpectrum::copy_members(const GCTAModelSpatialGaussSpectrum& model) +{ + // Clone sigma spectrum + if (model.m_sigma != NULL) { + m_sigma = model.m_sigma->clone(); + } + + // Set parameter pointers + set_pointers(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void GCTAModelSpatialGaussSpectrum::free_members(void) +{ + // Delete sigma spectrum + if (m_sigma != NULL) delete m_sigma; + + // Signal free pointer + m_sigma = NULL; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set pointers + * + * Set pointers to all model parameters. The pointers are stored in a vector + * that is member of the GModelData base class. + ***************************************************************************/ +void GCTAModelSpatialGaussSpectrum::set_pointers(void) +{ + // Clear parameters + m_pars.clear(); + + // Determine number of sigma parameters + int n_sigma = (m_sigma != NULL) ? m_sigma->size() : 0; + + // Gather sigma parameters if there are some + if (n_sigma > 0) { + for (int i = 0; i < n_sigma; ++i) { + m_pars.push_back(&((*m_sigma)[i])); + } + } + + // Return + return; +} diff --git a/inst/cta/test/data/cta_model_bgd_gauss_e.xml b/inst/cta/test/data/cta_model_bgd_gauss_e.xml new file mode 100644 index 0000000..6588998 --- /dev/null +++ b/inst/cta/test/data/cta_model_bgd_gauss_e.xml @@ -0,0 +1,22 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index d7ee04c..8170940 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -59,6 +59,7 @@ const std::string cta_unbin_xml = datadir+"/obs_unbinned.xml"; const std::string cta_model_xml = datadir+"/crab.xml"; const std::string cta_rsp_xml = datadir+"/rsp_models.xml"; const std::string cta_bgd_rad_gauss_xml = datadir+"/cta_model_bgd_rad_gauss.xml"; +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"; @@ -176,6 +177,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_gauss_spectrum), + "Test GCTAModelSpatialGaussSpectrum class"); // Return return; @@ -2029,6 +2032,57 @@ void TestGCTAModel::test_model_aeff_bgd(void) /***********************************************************************//** + * @brief Test GCTAModelSpatialGaussSpectrum class + ***************************************************************************/ +void TestGCTAModel::test_spatial_gauss_spectrum(void) +{ + // Set instrument direction, energy and time + GCTAInstDir dir(0.0, 0.0); + GEnergy energy(1.0, "TeV"); + GTime time; + + // Test void constructor + GCTAModelSpatialGaussSpectrum model1; + test_assert(model1.sigma() == NULL, "Test sigma() for void constructor"); + test_value(model1.eval(dir, energy, time), 0.0, "Test eval() for void constructor"); + + // Test value constructor + GCTAModelSpatialGaussSpectrum model2(3.0); + test_assert(model2.sigma() != NULL, "Test sigma() for value constructor"); + test_value(model2.eval(dir, energy, time), 1.0, "Test eval() for value constructor"); + + // Test spectrum constructor + GModelSpectralConst spectrum(3.0); + GCTAModelSpatialGaussSpectrum model3(spectrum); + test_assert(model3.sigma() != NULL, "Test sigma() for spectrum constructor"); + test_value(model3.eval(dir, energy, time), 1.0, "Test eval() for spectrum constructor"); + + // Test XML constructor + GXml xml(cta_bgd_gauss_e_xml); + GXmlElement* element = xml.element("source_library > source[0] > spatialModel"); + GCTAModelSpatialGaussSpectrum model4(*element); + test_assert(model4.sigma() != NULL, "Test sigma() for XML constructor"); + test_value(model4.eval(dir, energy, time), 1.0, "Test eval() for XML constructor"); + + // Test copy constructor + GCTAModelSpatialGaussSpectrum model5(model4); + test_assert(model5.sigma() != NULL, "Test sigma() for copy constructor"); + test_value(model5.eval(dir, energy, time), 1.0, "Test eval() for copy constructor"); + + // Test sigma() setter + GCTAModelSpatialGaussSpectrum model6; + test_assert(model6.sigma() == NULL, "Test sigma() setter"); + test_value(model6.eval(dir, energy, time), 0.0, "Test sigma() setter"); + model6.sigma(spectrum); + test_assert(model6.sigma() != NULL, "Test sigma() setter"); + test_value(model6.eval(dir, energy, time), 1.0, "Test sigma() setter"); + + // Return + return; +} + + +/***********************************************************************//** * @brief Checks handling of CTA event bin ***************************************************************************/ void TestGCTAObservation::test_event_bin(void) diff --git a/inst/cta/test/test_CTA.hpp b/inst/cta/test/test_CTA.hpp index 6272d82..2e6cd5c 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_gauss_spectrum(void); };