Olga Sergijenko / gammalib

Commit 4f9fa41aaa185ed8c595a2a59cb85b4f53c332f6

Authored by Jürgen Knödlseder
1 parent 20f137dc

Add GCTAModelSkyCube class (#3442)

Inline Side-by-side
Showing 19 changed files with 1561 additions and 42 deletions
ChangeLog
  View file @4f9fa41
1   -2020-11-08
  1 +2020-11-09
2 2  
3 3 * Version 2.0.0 released
4 4 ========================
5 5  
  6 + Add GCTAModelSkyCube class (#3442)
6 7 Add GSkyModel flux() and eflux() methods (#3439)
7 8 Add XML value constructors and access methods to GXmlElement class (#3431)
8 9 Correct application elapsed CPU time for OpenMP (#3417)
... ...
1 1 New Features and Important Changes in GammaLib 2.0.0
2 2  
3   -8 November 2020
  3 +9 November 2020
4 4  
5 5  
6 6 1. Introduction
... ... @@ -15,6 +15,7 @@ The following classes have been added:
15 15 - GFunctions
16 16 - GIntegrals
17 17 - GResponseVectorCache
  18 +- GCTAModelSkyCube
18 19  
19 20 The following classes have been removed:
20 21 - none
... ... @@ -181,6 +182,8 @@ Adapted the GLATResponse interface to the vectorised response computation (#3203
181 182  
182 183 20. CTA interface
183 184 -----------------
  185 +Added the GCTAModelSkyCube class that implements a data space model for CTA (#3442).
  186 +
184 187 Implement in GCTAResponseCube an optimised vectorised response computation for radial
185 188 models that applies for a stacked analysis. The code supports in particular analytical
186 189 spatial model gradients (#3203).
... ...
README.md
  View file @4f9fa41
1 1 GammaLib information
2 2 ====================
3   -* Version: 2.0.0.dev (8 November 2020)
  3 +* Version: 2.0.0.dev (9 November 2020)
4 4  
5 5 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/)
6 6  
... ...
doc/source/admin/release_history/2.0.rst
  View file @4f9fa41
1 1 .. _2.0:
2 2  
3   -GammaLib 2.0 release (8 November 2020)
  3 +GammaLib 2.0 release (9 November 2020)
4 4 ======================================
5 5  
6 6 Overview
... ... @@ -24,18 +24,12 @@ Bug fixes
24 24 Improvements
25 25 ------------
26 26  
27   -* [`3439 <https://cta-redmine.irap.omp.eu/issues/3439>`_] -
28   - Add ``GSkyModel`` ``flux()`` and ``eflux()`` methods
29   -* [`3431 <https://cta-redmine.irap.omp.eu/issues/3431>`_] -
30   - Add XML value constructors and access methods to ``GXmlElement`` class
31 27 * [`3333 <https://cta-redmine.irap.omp.eu/issues/3333>`_] -
32 28 Integrate analytical gradients for energy dispersion computation
33 29 * Implement ``GMatrixSparse::row()`` and ``GMatrixSparse::add_to_row()`` methods
34 30 * [`3203 <https://cta-redmine.irap.omp.eu/issues/3203>`_] -
35 31 Implement vectorised response computation and analytical gradients for
36 32 the radial Gaussian model in a stacked analysis
37   -* [`3074 <https://cta-redmine.irap.omp.eu/issues/3074>`_] -
38   - Add ``GSkyModel`` ``flux_error()`` and ``eflux_error()`` methods
39 33 * [`1300 <https://cta-redmine.irap.omp.eu/issues/1300>`_] -
40 34 Implement pre-computation cache in GModelSpatialGauss
41 35  
... ... @@ -43,10 +37,18 @@ Improvements
43 37 New features
44 38 ------------
45 39  
  40 +* [`3442 <https://cta-redmine.irap.omp.eu/issues/3442>`_] -
  41 + Add ``GCTAModelSkyCube`` class
  42 +* [`3439 <https://cta-redmine.irap.omp.eu/issues/3439>`_] -
  43 + Add ``GSkyModel`` ``flux()`` and ``eflux()`` methods
  44 +* [`3431 <https://cta-redmine.irap.omp.eu/issues/3431>`_] -
  45 + Add XML value constructors and access methods to ``GXmlElement`` class
46 46 * [`3403 <https://cta-redmine.irap.omp.eu/issues/3403>`_] -
47 47 Implement sparse vector response cache
48 48 * [`3328 <https://cta-redmine.irap.omp.eu/issues/3328>`_] -
49 49 Add ``GFunctions`` abstract base class and ``GIntegrals`` class
  50 +* [`3074 <https://cta-redmine.irap.omp.eu/issues/3074>`_] -
  51 + Add ``GSkyModel`` ``flux_error()`` and ``eflux_error()`` methods
50 52  
51 53  
52 54 Clean up
... ...
inst/cta/Makefile.am
  View file @4f9fa41
... ... @@ -88,7 +88,8 @@ sources = src/cta_helpers_response_stacked_vector.cpp \
88 88 src/GCTAModelRadialGauss.cpp \
89 89 src/GCTAModelRadialPolynom.cpp \
90 90 src/GCTAModelRadialProfile.cpp \
91   - src/GCTAModelRadialAcceptance.cpp
  91 + src/GCTAModelRadialAcceptance.cpp \
  92 + src/GCTAModelSkyCube.cpp
92 93  
93 94 # Define headers to be installed
94 95 pkginclude_HEADERS = include/GCTATypemaps.hpp \
... ... @@ -144,6 +145,7 @@ pkginclude_HEADERS = include/GCTATypemaps.hpp \
144 145 include/GCTAModelRadialPolynom.hpp \
145 146 include/GCTAModelRadialProfile.hpp \
146 147 include/GCTAModelRadialAcceptance.hpp \
  148 + include/GCTAModelSkyCube.hpp \
147 149 include/GCTALib.hpp
148 150  
149 151 # Build libtool library
... ... @@ -178,7 +180,8 @@ testdatacta = test/data/irf_1dc.xml \
178 180 test/data/crab_irf.xml \
179 181 test/data/onoff_model.xml \
180 182 test/data/crab_events.fits.gz \
181   - test/data/crab_cntmap.fits.gz
  183 + test/data/crab_cntmap.fits.gz \
  184 + test/data/crab_modcube.fits.gz
182 185 testdatactadir = $(pyexecdir)/gammalib/tests/cta/data
183 186 testdatacta_DATA = $(testdatacta)
184 187 testcaldbcta = test/caldb/cta_dummy_irf.dat \
... ...
inst/cta/include/GCTALib.hpp
  View file @4f9fa41
... ... @@ -86,5 +86,6 @@
86 86 #include "GCTAModelRadialPolynom.hpp"
87 87 #include "GCTAModelRadialProfile.hpp"
88 88 #include "GCTAModelRadialAcceptance.hpp"
  89 +#include "GCTAModelSkyCube.hpp"
89 90  
90 91 #endif /* GCTALIB_HPP */
... ...
inst/cta/include/GCTAModelSkyCube.hpp 0 → 100644
  View file @4f9fa41
  1 +/***************************************************************************
  2 + * GCTAModelSkyCube.hpp - CTA sky cube model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2020 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 GCTAModelSkyCube.hpp
  23 + * @brief CTA sky cube model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GCTAMODELSKYCUBE_HPP
  28 +#define GCTAMODELSKYCUBE_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include <vector>
  33 +#include "GModelData.hpp"
  34 +#include "GFilename.hpp"
  35 +#include "GSkyMap.hpp"
  36 +#include "GEbounds.hpp"
  37 +#include "GNodeArray.hpp"
  38 +#include "GModelPar.hpp"
  39 +#include "GModelSpectral.hpp"
  40 +#include "GModelTemporal.hpp"
  41 +#include "GCTAEventList.hpp"
  42 +
  43 +/* __ Forward declarations _______________________________________________ */
  44 +class GEvent;
  45 +class GObservation;
  46 +class GXmlElement;
  47 +
  48 +
  49 +/***********************************************************************//**
  50 + * @class GCTAModelSkyCube
  51 + *
  52 + * @brief CTA sky cube model class
  53 + *
  54 + * This class implements a sky cube model for CTA.
  55 + ***************************************************************************/
  56 +class GCTAModelSkyCube : public GModelData {
  57 +
  58 +public:
  59 + // Constructors and destructors
  60 + GCTAModelSkyCube(void);
  61 + explicit GCTAModelSkyCube(const GXmlElement& xml);
  62 + GCTAModelSkyCube(const GFilename& filename,
  63 + const GModelSpectral& spectral);
  64 + GCTAModelSkyCube(const GFilename& filename,
  65 + const GModelSpectral& spectral,
  66 + const GModelTemporal& temporal);
  67 + GCTAModelSkyCube(const GCTAModelSkyCube& bgd);
  68 + virtual ~GCTAModelSkyCube(void);
  69 +
  70 + // Operators
  71 + virtual GCTAModelSkyCube& operator=(const GCTAModelSkyCube& model);
  72 +
  73 + // Implemented pure virtual methods
  74 + virtual void clear(void);
  75 + virtual GCTAModelSkyCube* clone(void) const;
  76 + virtual std::string classname(void) const;
  77 + virtual std::string type(void) const;
  78 + virtual bool is_constant(void) const;
  79 + virtual double eval(const GEvent& event,
  80 + const GObservation& obs,
  81 + const bool& gradients = false) const;
  82 + virtual double npred(const GEnergy& obsEng,
  83 + const GTime& obsTime,
  84 + const GObservation& obs) const;
  85 + virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
  86 + virtual void read(const GXmlElement& xml);
  87 + virtual void write(GXmlElement& xml) const;
  88 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  89 +
  90 + // Other methods
  91 + void load(const GFilename& filename);
  92 + const GFilename& filename(void) const;
  93 + GModelSpectral* spectral(void) const;
  94 + GModelTemporal* temporal(void) const;
  95 + void spectral(const GModelSpectral* spectral);
  96 + void temporal(const GModelTemporal* temporal);
  97 +
  98 +protected:
  99 + // Methods
  100 + void init_members(void);
  101 + void copy_members(const GCTAModelSkyCube& bgd);
  102 + void free_members(void);
  103 + void set_pointers(void);
  104 + bool valid_model(void) const;
  105 + void read_xml_spatial(const GXmlElement& xml);
  106 + void write_xml_spatial(GXmlElement& xml) const;
  107 + GModelSpectral* xml_spectral(const GXmlElement& spectral) const;
  108 + GModelTemporal* xml_temporal(const GXmlElement& temporal) const;
  109 +
  110 + // Members
  111 + GFilename m_filename; //!< Filename
  112 + GSkyMap m_cube; //!< Sky cube
  113 + GEbounds m_ebounds; //!< Energy boundaries of sky cube
  114 + GNodeArray m_elogmeans; //!< Mean log10(TeV) energies for the sky cube
  115 + GModelPar m_norm; //!< Normalisation parameter
  116 + GModelSpectral* m_spectral; //!< Spectral model
  117 + GModelTemporal* m_temporal; //!< Temporal model
  118 +
  119 + // Npred cache
  120 + mutable std::vector<std::string> m_npred_names; //!< Model names
  121 + mutable std::vector<GEnergy> m_npred_energies; //!< Model energy
  122 + mutable std::vector<GTime> m_npred_times; //!< Model time
  123 + mutable std::vector<double> m_npred_values; //!< Model values
  124 +};
  125 +
  126 +
  127 +/***********************************************************************//**
  128 + * @brief Return class name
  129 + *
  130 + * @return String containing the class name ("GCTAModelSkyCube").
  131 + ***************************************************************************/
  132 +inline
  133 +std::string GCTAModelSkyCube::classname(void) const
  134 +{
  135 + return ("GCTAModelSkyCube");
  136 +}
  137 +
  138 +
  139 +/***********************************************************************//**
  140 + * @brief Return data model type
  141 + *
  142 + * @return Data model type.
  143 + *
  144 + * Returns the type of the data model.
  145 + ***************************************************************************/
  146 +inline
  147 +std::string GCTAModelSkyCube::type(void) const
  148 +{
  149 + return ("CTASkyCube");
  150 +}
  151 +
  152 +
  153 +/***********************************************************************//**
  154 + * @brief Signals if sky model is temporally constant
  155 + *
  156 + * @return True if sky model is temporally constant, false otherwise.
  157 + *
  158 + * Signals if the sky model is temporally constant. A temporally constant
  159 + * model is a model that has a temporal component of type "Constant". If
  160 + * the model has no temporal component it is also consider being constant.
  161 + ***************************************************************************/
  162 +inline
  163 +bool GCTAModelSkyCube::is_constant(void) const
  164 +{
  165 + return ((m_temporal != NULL && m_temporal->type() == "Constant") ||
  166 + (m_temporal == NULL));
  167 +}
  168 +
  169 +
  170 +/***********************************************************************//**
  171 + * @brief Return filename
  172 + *
  173 + * @return Filename.
  174 + ***************************************************************************/
  175 +inline
  176 +const GFilename& GCTAModelSkyCube::filename(void) const
  177 +{
  178 + return (m_filename);
  179 +}
  180 +
  181 +
  182 +/***********************************************************************//**
  183 + * @brief Return spectral model component
  184 + *
  185 + * @return Pointer to spectral model component.
  186 + *
  187 + * Returns a pointer to the spectral model component of the model. The
  188 + * pointer is of type GModelSpectral. Note that a NULL pointer may be
  189 + * returned if the sky model has no spectral model component.
  190 + ***************************************************************************/
  191 +inline
  192 +GModelSpectral* GCTAModelSkyCube::spectral(void) const
  193 +{
  194 + return (m_spectral);
  195 +}
  196 +
  197 +
  198 +/***********************************************************************//**
  199 + * @brief Return temporal model component
  200 + *
  201 + * @return Pointer to temporal model component.
  202 + *
  203 + * Returns a pointer to the temporal model component of the model. The
  204 + * pointer is of type GModelTemporal. Note that a NULL pointer may be
  205 + * returned if the sky model has no temporal model component.
  206 + ***************************************************************************/
  207 +inline
  208 +GModelTemporal* GCTAModelSkyCube::temporal(void) const
  209 +{
  210 + return (m_temporal);
  211 +}
  212 +
  213 +#endif /* GCTAMODELSKYCUBE_HPP */
... ...
inst/cta/pyext/GCTAModelSkyCube.i 0 → 100644
  View file @4f9fa41
  1 +/***************************************************************************
  2 + * GCTAModelSkyCube.i - CTA sky cube model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2020 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 GCTAModelSkyCube.i
  23 + * @brief CTA sky cube model class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCTAModelSkyCube.hpp"
  29 +#include "GTools.hpp"
  30 +%}
  31 +
  32 +
  33 +/***********************************************************************//**
  34 +* @class GCTAModelSkyCube
  35 +*
  36 +* @brief CTA sky cube model class
  37 +***************************************************************************/
  38 +class GCTAModelSkyCube : public GModelData {
  39 +public:
  40 + // Constructors and destructors
  41 + GCTAModelSkyCube(void);
  42 + explicit GCTAModelSkyCube(const GXmlElement& xml);
  43 + GCTAModelSkyCube(const GFilename& filename,
  44 + const GModelSpectral& spectral);
  45 + GCTAModelSkyCube(const GFilename& filename,
  46 + const GModelSpectral& spectral,
  47 + const GModelTemporal& temporal);
  48 + GCTAModelSkyCube(const GCTAModelSkyCube& bgd);
  49 + virtual ~GCTAModelSkyCube(void);
  50 +
  51 + // Implemented pure virtual methods
  52 + virtual void clear(void);
  53 + virtual GCTAModelSkyCube* clone(void) const;
  54 + virtual std::string classname(void) const;
  55 + virtual std::string type(void) const;
  56 + virtual bool is_constant(void) const;
  57 + virtual double eval(const GEvent& event,
  58 + const GObservation& obs,
  59 + const bool& gradients = false) const;
  60 + virtual double npred(const GEnergy& obsEng,
  61 + const GTime& obsTime,
  62 + const GObservation& obs) const;
  63 + virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
  64 + virtual void read(const GXmlElement& xml);
  65 + virtual void write(GXmlElement& xml) const;
  66 +
  67 + // Other methods
  68 + void load(const GFilename& filename);
  69 + const GFilename& filename(void) const;
  70 + GModelSpectral* spectral(void) const;
  71 + GModelTemporal* temporal(void) const;
  72 + void spectral(const GModelSpectral* spectral);
  73 + void temporal(const GModelTemporal* temporal);
  74 +};
  75 +
  76 +/***********************************************************************//**
  77 +* @brief GCTAModelSkyCube class extension
  78 +***************************************************************************/
  79 +%extend GCTAModelSkyCube {
  80 + GCTAModelSkyCube copy() {
  81 + return (*self);
  82 + }
  83 +%pythoncode {
  84 + def __getstate__(self):
  85 + state = (self.filename(), self.spectral(), self.temporal(),
  86 + gammalib.GModelData.__getstate__(self))
  87 + return state
  88 + def __setstate__(self, state):
  89 + if state[0] == '' and state[1] == None and state[2] == None:
  90 + self.__init__()
  91 + elif state[2] == None:
  92 + self.__init__(state[0], state[1])
  93 + else:
  94 + self.__init__(state[0], state[1], state[2])
  95 + gammalib.GModelData.__setstate__(self, state[3])
  96 +}
  97 +};
... ...
inst/cta/pyext/cta.i
  View file @4f9fa41
... ... @@ -403,3 +403,4 @@
403 403 %include "GCTAModelRadialPolynom.i"
404 404 %include "GCTAModelRadialProfile.i"
405 405 %include "GCTAModelRadialAcceptance.i"
  406 +%include "GCTAModelSkyCube.i"
... ...
inst/cta/src/GCTAModelAeffBackground.cpp
  View file @4f9fa41
... ... @@ -754,9 +754,9 @@ GCTAEventList* GCTAModelAeffBackground::mc(const GObservation&amp; obs,
754 754 * <spectrum type="...">
755 755 * ...
756 756 * </spectrum>
757   - * <temporalModel type="...">
  757 + * <temporal type="...">
758 758 * ...
759   - * </temporalModel>
  759 + * </temporal>
760 760 * </source>
761 761 *
762 762 * If no temporal component is found a constant model is assumed.
... ... @@ -776,8 +776,8 @@ void GCTAModelAeffBackground::read(const GXmlElement&amp; xml)
776 776 m_spectral = xml_spectral(*spectral);
777 777  
778 778 // Optionally get temporal model
779   - if (xml.elements("temporalModel")) {
780   - const GXmlElement* temporal = xml.element("temporalModel", 0);
  779 + if (xml.elements("temporal")) {
  780 + const GXmlElement* temporal = xml.element("temporal", 0);
781 781 m_temporal = xml_temporal(*temporal);
782 782 }
783 783 else {
... ... @@ -817,9 +817,9 @@ void GCTAModelAeffBackground::read(const GXmlElement&amp; xml)
817 817 * <spectrum type="...">
818 818 * ...
819 819 * </spectrum>
820   - * <temporalModel type="...">
  820 + * <temporal type="...">
821 821 * ...
822   - * </temporalModel>
  822 + * </temporal>
823 823 * </source>
824 824 *
825 825 * If no temporal component is found a constant model is assumed.
... ... @@ -852,7 +852,7 @@ void GCTAModelAeffBackground::write(GXmlElement&amp; xml) const
852 852 if (src == NULL) {
853 853 src = xml.append("source");
854 854 if (spectral() != NULL) src->append(GXmlElement("spectrum"));
855   - if (write_temporal) src->append(GXmlElement("temporalModel"));
  855 + if (write_temporal) src->append(GXmlElement("temporal"));
856 856 }
857 857  
858 858 // Write spectral model
... ... @@ -864,7 +864,7 @@ void GCTAModelAeffBackground::write(GXmlElement&amp; xml) const
864 864 // Optionally write temporal model
865 865 if (write_temporal) {
866 866 if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
867   - GXmlElement* temp = src->element("temporalModel", 0);
  867 + GXmlElement* temp = src->element("temporal", 0);
868 868 temporal()->write(*temp);
869 869 }
870 870 }
... ...
inst/cta/src/GCTAModelCubeBackground.cpp
  View file @4f9fa41
... ... @@ -537,9 +537,9 @@ GCTAEventList* GCTAModelCubeBackground::mc(const GObservation&amp; obs, GRan&amp; ran) c
537 537 * <spectrum type="...">
538 538 * ...
539 539 * </spectrum>
540   - * <temporalModel type="...">
  540 + * <temporal type="...">
541 541 * ...
542   - * </temporalModel>
  542 + * </temporal>
543 543 * </source>
544 544 *
545 545 * If no temporal component is found a constant model is assumed.
... ... @@ -559,8 +559,8 @@ void GCTAModelCubeBackground::read(const GXmlElement&amp; xml)
559 559 m_spectral = xml_spectral(*spectral);
560 560  
561 561 // Optionally get temporal model
562   - if (xml.elements("temporalModel")) {
563   - const GXmlElement* temporal = xml.element("temporalModel", 0);
  562 + if (xml.elements("temporal")) {
  563 + const GXmlElement* temporal = xml.element("temporal", 0);
564 564 m_temporal = xml_temporal(*temporal);
565 565 }
566 566 else {
... ... @@ -600,9 +600,9 @@ void GCTAModelCubeBackground::read(const GXmlElement&amp; xml)
600 600 * <spectrum type="...">
601 601 * ...
602 602 * </spectrum>
603   - * <temporalModel type="...">
  603 + * <temporal type="...">
604 604 * ...
605   - * </temporalModel>
  605 + * </temporal>
606 606 * </source>
607 607 *
608 608 * If no temporal component is found a constant model is assumed.
... ... @@ -635,7 +635,7 @@ void GCTAModelCubeBackground::write(GXmlElement&amp; xml) const
635 635 if (src == NULL) {
636 636 src = xml.append("source");
637 637 if (spectral() != NULL) src->append(GXmlElement("spectrum"));
638   - if (write_temporal) src->append(GXmlElement("temporalModel"));
  638 + if (write_temporal) src->append(GXmlElement("temporal"));
639 639 }
640 640  
641 641 // Write spectral model
... ... @@ -647,7 +647,7 @@ void GCTAModelCubeBackground::write(GXmlElement&amp; xml) const
647 647 // Optionally write temporal model
648 648 if (write_temporal) {
649 649 if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
650   - GXmlElement* temp = src->element("temporalModel", 0);
  650 + GXmlElement* temp = src->element("temporal", 0);
651 651 temporal()->write(*temp);
652 652 }
653 653 }
... ...
inst/cta/src/GCTAModelIrfBackground.cpp
  View file @4f9fa41
... ... @@ -685,9 +685,9 @@ GCTAEventList* GCTAModelIrfBackground::mc(const GObservation&amp; obs, GRan&amp; ran) co
685 685 * <spectrum type="...">
686 686 * ...
687 687 * </spectrum>
688   - * <temporalModel type="...">
  688 + * <temporal type="...">
689 689 * ...
690   - * </temporalModel>
  690 + * </temporal>
691 691 * </source>
692 692 *
693 693 * If no temporal component is found a constant model is assumed.
... ... @@ -707,8 +707,8 @@ void GCTAModelIrfBackground::read(const GXmlElement&amp; xml)
707 707 m_spectral = xml_spectral(*spectral);
708 708  
709 709 // Optionally get temporal model
710   - if (xml.elements("temporalModel")) {
711   - const GXmlElement* temporal = xml.element("temporalModel", 0);
  710 + if (xml.elements("temporal")) {
  711 + const GXmlElement* temporal = xml.element("temporal", 0);
712 712 m_temporal = xml_temporal(*temporal);
713 713 }
714 714 else {
... ... @@ -748,9 +748,9 @@ void GCTAModelIrfBackground::read(const GXmlElement&amp; xml)
748 748 * <spectrum type="...">
749 749 * ...
750 750 * </spectrum>
751   - * <temporalModel type="...">
  751 + * <temporal type="...">
752 752 * ...
753   - * </temporalModel>
  753 + * </temporal>
754 754 * </source>
755 755 *
756 756 * If no temporal component is found a constant model is assumed.
... ... @@ -783,7 +783,7 @@ void GCTAModelIrfBackground::write(GXmlElement&amp; xml) const
783 783 if (src == NULL) {
784 784 src = xml.append("source");
785 785 if (spectral() != NULL) src->append(GXmlElement("spectrum"));
786   - if (write_temporal) src->append(GXmlElement("temporalModel"));
  786 + if (write_temporal) src->append(GXmlElement("temporal"));
787 787 }
788 788  
789 789 // Write spectral model
... ... @@ -795,7 +795,7 @@ void GCTAModelIrfBackground::write(GXmlElement&amp; xml) const
795 795 // Optionally write temporal model
796 796 if (write_temporal) {
797 797 if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
798   - GXmlElement* temp = src->element("temporalModel", 0);
  798 + GXmlElement* temp = src->element("temporal", 0);
799 799 temporal()->write(*temp);
800 800 }
801 801 }
... ...
inst/cta/src/GCTAModelRadialAcceptance.cpp
  View file @4f9fa41
... ... @@ -653,7 +653,7 @@ GCTAEventList* GCTAModelRadialAcceptance::mc(const GObservation&amp; obs,
653 653 *
654 654 * The model is composed of a spectrum component ('spectral'), a radial
655 655 * component ('radialModel'), and, optionally, of a temporal component
656   - * ('lightcurve'). If no temporal component is found a constant model is
  656 + * ('temporal'). If no temporal component is found a constant model is
657 657 * assumed.
658 658 ***************************************************************************/
659 659 void GCTAModelRadialAcceptance::read(const GXmlElement& xml)
... ... @@ -675,7 +675,7 @@ void GCTAModelRadialAcceptance::read(const GXmlElement&amp; xml)
675 675  
676 676 // Optionally get temporal model
677 677 try {
678   - const GXmlElement* temp = xml.element("lightcurve", 0);
  678 + const GXmlElement* temp = xml.element("temporal", 0);
679 679 m_temporal = xml_temporal(*temp);
680 680 }
681 681 catch (GException::xml_name_not_found &e) {
... ... @@ -730,7 +730,7 @@ void GCTAModelRadialAcceptance::write(GXmlElement&amp; xml) const
730 730 src = xml.append("source");
731 731 if (spectral() != NULL) src->append(GXmlElement("spectrum"));
732 732 if (radial() != NULL) src->append(GXmlElement("radialModel"));
733   - if (write_temporal) src->append(GXmlElement("temporalModel"));
  733 + if (write_temporal) src->append(GXmlElement("temporal"));
734 734 }
735 735  
736 736 // Write radial model
... ... @@ -748,7 +748,7 @@ void GCTAModelRadialAcceptance::write(GXmlElement&amp; xml) const
748 748 // Optionally write temporal model
749 749 if (write_temporal) {
750 750 if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
751   - GXmlElement* temp = src->element("temporalModel", 0);
  751 + GXmlElement* temp = src->element("temporal", 0);
752 752 temporal()->write(*temp);
753 753 }
754 754 }
... ...
inst/cta/src/GCTAModelSkyCube.cpp 0 → 100644
  View file @4f9fa41
  1 +/***************************************************************************
  2 + * GCTAModelSkyCube.cpp - CTA sky cube model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2020 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 GCTAModelSkyCube.cpp
  23 + * @brief CTA sky cube model class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GModelRegistry.hpp"
  33 +#include "GModelSpectralRegistry.hpp"
  34 +#include "GModelTemporalRegistry.hpp"
  35 +#include "GModelTemporalConst.hpp"
  36 +#include "GObservation.hpp"
  37 +#include "GCTAModelSkyCube.hpp"
  38 +#include "GCTASupport.hpp"
  39 +
  40 +/* __ Globals ____________________________________________________________ */
  41 +const GCTAModelSkyCube g_cta_inst_sky_seed;
  42 +const GModelRegistry g_cta_inst_sky_registry(&g_cta_inst_sky_seed);
  43 +
  44 +/* __ Method name definitions ____________________________________________ */
  45 +#define G_EVAL "GCTAModelSkyCube::eval(GEvent&, GObservation&, bool&)"
  46 +#define G_NPRED "GCTAModelSkyCube::npred(GEnergy&, GTime&, GObservation&)"
  47 +#define G_MC "GCTAModelSkyCube::mc(GObservation&, GRan&)"
  48 +#define G_READ_XML_SPATIAL "GCTAModelSkyCube::read_xml_spatial(GXmlElement&)"
  49 +#define G_WRITE_XML_SPATIAL "GCTAModelSkyCube::write_xml_spatial("\
  50 + "GXmlElement&)"
  51 +#define G_XML_SPECTRAL "GCTAModelSkyCube::xml_spectral(GXmlElement&)"
  52 +#define G_XML_TEMPORAL "GCTAModelSkyCube::xml_temporal(GXmlElement&)"
  53 +
  54 +/* __ Macros _____________________________________________________________ */
  55 +
  56 +/* __ Coding definitions _________________________________________________ */
  57 +#define G_USE_NPRED_CACHE
  58 +
  59 +/* __ Debug definitions __________________________________________________ */
  60 +//#define G_DEBUG_NPRED
  61 +
  62 +/* __ Constants __________________________________________________________ */
  63 +
  64 +
  65 +/*==========================================================================
  66 + = =
  67 + = Constructors/destructors =
  68 + = =
  69 + ==========================================================================*/
  70 +
  71 +/***********************************************************************//**
  72 + * @brief Void constructor
  73 + ***************************************************************************/
  74 +GCTAModelSkyCube::GCTAModelSkyCube(void) : GModelData()
  75 +{
  76 + // Initialise class members
  77 + init_members();
  78 +
  79 + // Return
  80 + return;
  81 +}
  82 +
  83 +
  84 +/***********************************************************************//**
  85 + * @brief XML constructor
  86 + *
  87 + * @param[in] xml XML element.
  88 + *
  89 + * Constructs a CTA sky cube model from the information provided by an
  90 + * XML element (see GCTAModelSkyCube::read method).
  91 + ***************************************************************************/
  92 +GCTAModelSkyCube::GCTAModelSkyCube(const GXmlElement& xml) : GModelData(xml)
  93 +{
  94 + // Initialise members
  95 + init_members();
  96 +
  97 + // Read XML
  98 + read(xml);
  99 +
  100 + // Set parameter pointers
  101 + set_pointers();
  102 +
  103 + // Return
  104 + return;
  105 +}
  106 +
  107 +
  108 +/***********************************************************************//**
  109 + * @brief Construct from filename and spectral component
  110 + *
  111 + * @param[in] filename Cube file name.
  112 + * @param[in] spectral Spectral model component.
  113 + *
  114 + * Constructs a CTA sky cube model from a cube @p filename and a @p spectral
  115 + * model component. The temporal component is assumed to be constant.
  116 + ***************************************************************************/
  117 +GCTAModelSkyCube::GCTAModelSkyCube(const GFilename& filename,
  118 + const GModelSpectral& spectral) :
  119 + GModelData()
  120 +{
  121 + // Initialise members
  122 + init_members();
  123 +
  124 + // Allocate temporal constant model
  125 + GModelTemporalConst temporal;
  126 +
  127 + // Load filename
  128 + load(filename);
  129 +
  130 + // Clone model components
  131 + m_spectral = spectral.clone();
  132 + m_temporal = temporal.clone();
  133 +
  134 + // Set parameter pointers
  135 + set_pointers();
  136 +
  137 + // Return
  138 + return;
  139 +}
  140 +
  141 +
  142 +/***********************************************************************//**
  143 + * @brief Construct from model components
  144 + *
  145 + * @param[in] filename Cube file name.
  146 + * @param[in] spectral Spectral model component.
  147 + * @param[in] temporal Temporal model component.
  148 + *
  149 + * Constructs a CTA sky cube model from a cube @p filename, a @p spectral
  150 + * model component and a @p temporal component.
  151 + ***************************************************************************/
  152 +GCTAModelSkyCube::GCTAModelSkyCube(const GFilename& filename,
  153 + const GModelSpectral& spectral,
  154 + const GModelTemporal& temporal) :
  155 + GModelData()
  156 +{
  157 + // Initialise members
  158 + init_members();
  159 +
  160 + // Load filename
  161 + load(filename);
  162 +
  163 + // Clone model components
  164 + m_spectral = spectral.clone();
  165 + m_temporal = temporal.clone();
  166 +
  167 + // Set parameter pointers
  168 + set_pointers();
  169 +
  170 + // Return
  171 + return;
  172 +}
  173 +
  174 +
  175 +/***********************************************************************//**
  176 + * @brief Copy constructor
  177 + *
  178 + * @param[in] sky CTA sky cube model.
  179 + ***************************************************************************/
  180 +GCTAModelSkyCube::GCTAModelSkyCube(const GCTAModelSkyCube& sky) :
  181 + GModelData(sky)
  182 +{
  183 + // Initialise class members
  184 + init_members();
  185 +
  186 + // Copy members
  187 + copy_members(sky);
  188 +
  189 + // Return
  190 + return;
  191 +}
  192 +
  193 +
  194 +/***********************************************************************//**
  195 + * @brief Destructor
  196 + ***************************************************************************/
  197 +GCTAModelSkyCube::~GCTAModelSkyCube(void)
  198 +{
  199 + // Free members
  200 + free_members();
  201 +
  202 + // Return
  203 + return;
  204 +}
  205 +
  206 +
  207 +/*==========================================================================
  208 + = =
  209 + = Operators =
  210 + = =
  211 + ==========================================================================*/
  212 +
  213 +/***********************************************************************//**
  214 + * @brief Assignment operator
  215 + *
  216 + * @param[in] sky CTA sky cube model.
  217 + * @return CTA sky cube model.
  218 + ***************************************************************************/
  219 +GCTAModelSkyCube& GCTAModelSkyCube::operator=(const GCTAModelSkyCube& sky)
  220 +{
  221 + // Execute only if object is not identical
  222 + if (this != &sky) {
  223 +
  224 + // Copy base class members
  225 + this->GModelData::operator=(sky);
  226 +
  227 + // Free members
  228 + free_members();
  229 +
  230 + // Initialise private members
  231 + init_members();
  232 +
  233 + // Copy members
  234 + copy_members(sky);
  235 +
  236 + } // endif: object was not identical
  237 +
  238 + // Return this object
  239 + return *this;
  240 +}
  241 +
  242 +
  243 +/*==========================================================================
  244 + = =
  245 + = Public methods =
  246 + = =
  247 + ==========================================================================*/
  248 +
  249 +/***********************************************************************//**
  250 + * @brief Clear CTA sky cube model
  251 + *
  252 + * This method properly resets the CTA sky cube model to an initial state.
  253 + ***************************************************************************/
  254 +void GCTAModelSkyCube::clear(void)
  255 +{
  256 + // Free class members (base and derived classes, derived class first)
  257 + free_members();
  258 + this->GModelData::free_members();
  259 +
  260 + // Initialise members
  261 + this->GModelData::init_members();
  262 + init_members();
  263 +
  264 + // Return
  265 + return;
  266 +}
  267 +
  268 +
  269 +/***********************************************************************//**
  270 + * @brief Clone CTA sky cube model
  271 + *
  272 + * @return Pointer to deep copy of CTA sky cube model.
  273 + ***************************************************************************/
  274 +GCTAModelSkyCube* GCTAModelSkyCube::clone(void) const
  275 +{
  276 + return new GCTAModelSkyCube(*this);
  277 +}
  278 +
  279 +
  280 +/***********************************************************************//**
  281 + * @brief Load sky cube
  282 + *
  283 + * @param[in] filename Sky cube filename.
  284 + *
  285 + * Load the sky cube from the primary image extension.
  286 + ***************************************************************************/
  287 +void GCTAModelSkyCube::load(const GFilename& filename)
  288 +{
  289 + // Clear object
  290 + clear();
  291 +
  292 + // Put into OpenMP criticial zone
  293 + #pragma omp critical(GCTAModelSkyCube_load)
  294 + {
  295 + // Open FITS file
  296 + GFits fits(filename);
  297 +
  298 + // Get HDUs
  299 + const GFitsImage& hdu_skycube = *fits.image("Primary");
  300 + const GFitsTable& hdu_ebounds = *fits.table(gammalib::extname_ebounds);
  301 +
  302 + // Read sky cube
  303 + m_cube.read(hdu_skycube);
  304 +
  305 + // Read energy boundaries
  306 + m_ebounds.read(hdu_ebounds);
  307 +
  308 + // Set energy node array
  309 + for (int i = 0; i < m_ebounds.size(); ++i) {
  310 + m_elogmeans.append(m_ebounds.elogmean(i).log10TeV());
  311 + }
  312 +
  313 + // Close FITS file
  314 + fits.close();
  315 +
  316 + } // end of OpenMP critical zone
  317 +
  318 + // Store filename
  319 + m_filename = filename;
  320 +
  321 + // Return
  322 + return;
  323 +}
  324 +
  325 +
  326 +/***********************************************************************//**
  327 + * @brief Set spectral model component
  328 + *
  329 + * @param[in] spectral Pointer to spectral model component.
  330 + *
  331 + * Sets the spectral model component of the model.
  332 + ***************************************************************************/
  333 +void GCTAModelSkyCube::spectral(const GModelSpectral* spectral)
  334 +{
  335 + // Free spectral model component only if it differs from current
  336 + // component. This prevents unintential deallocation of the argument
  337 + if ((m_spectral != NULL) && (m_spectral != spectral)) {
  338 + delete m_spectral;
  339 + }
  340 +
  341 + // Clone spectral model component if it exists, otherwise set pointer
  342 + // to NULL
  343 + m_spectral = (spectral != NULL) ? spectral->clone() : NULL;
  344 +
  345 + // Set parameter pointers
  346 + set_pointers();
  347 +
  348 + // Return
  349 + return;
  350 +}
  351 +
  352 +
  353 +/***********************************************************************//**
  354 + * @brief Set temporal model component
  355 + *
  356 + * @param[in] temporal Pointer to temporal model component.
  357 + *
  358 + * Sets the temporal model component of the model.
  359 + ***************************************************************************/
  360 +void GCTAModelSkyCube::temporal(const GModelTemporal* temporal)
  361 +{
  362 + // Free temporal model component only if it differs from current
  363 + // component. This prevents unintential deallocation of the argument
  364 + if ((m_temporal != NULL) && (m_temporal != temporal)) {
  365 + delete m_temporal;
  366 + }
  367 +
  368 + // Clone temporal model component if it exists, otherwise set pointer
  369 + // to NULL
  370 + m_temporal = (temporal != NULL) ? temporal->clone() : NULL;
  371 +
  372 + // Set parameter pointers
  373 + set_pointers();
  374 +
  375 + // Return
  376 + return;
  377 +}
  378 +
  379 +
  380 +/***********************************************************************//**
  381 + * @brief Return event rate in units of events MeV\f$^{-1}\f$
  382 + * s\f$^{-1}\f$ sr\f$^{-1}\f$
  383 + *
  384 + * @param[in] event Observed event.
  385 + * @param[in] obs Observation.
  386 + * @param[in] gradients Compute gradients?
  387 + * @return Event rate (events MeV\f$^{-1}\f$ s\f$^{-1}\f$ sr\f$^{-1}\f$).
  388 + *
  389 + * Evaluates the sky cube model. The method returns a measured rate, defined
  390 + * as the number of counts per MeV, steradian and ontime.
  391 + *
  392 + * If the @p gradients flag is true the method will also set the parameter
  393 + * gradients of the model parameters.
  394 + ***************************************************************************/
  395 +double GCTAModelSkyCube::eval(const GEvent& event,
  396 + const GObservation& obs,
  397 + const bool& gradients) const
  398 +{
  399 + // Set indices and weighting factors for energy interpolation
  400 + m_elogmeans.set_value(event.energy().log10TeV());
  401 + int inx_left = m_elogmeans.inx_left();
  402 + int inx_right = m_elogmeans.inx_right();
  403 + double wgt_left = m_elogmeans.wgt_left();
  404 + double wgt_right = m_elogmeans.wgt_right();
  405 +
  406 + // Get reference on CTA instrument direction from event
  407 + const GCTAInstDir& dir = gammalib::cta_dir(G_EVAL, event);
  408 +
  409 + // Initialise non-normalised spatial component
  410 + double spat_raw = 0.0;
  411 +
  412 + // If left weight is close to 1, use left node
  413 + if (std::abs(wgt_left-1.0) < 1.0e-6) {
  414 + spat_raw = m_cube(dir.dir(), inx_left);
  415 + }
  416 +
  417 + // ... otherwise, if right weight is close to 1, use right node
  418 + else if (std::abs(wgt_right-1.0) < 1.0e-6) {
  419 + spat_raw = m_cube(dir.dir(), inx_right);
  420 + }
  421 +
  422 + // ... otherwise perform interpolation
  423 + else {
  424 + double spat_left = m_cube(dir.dir(), inx_left);
  425 + double spat_right = m_cube(dir.dir(), inx_right);
  426 + if (spat_left > 0.0 && spat_right > 0.0) {
  427 + spat_raw = std::exp(wgt_left * std::log(spat_left) +
  428 + wgt_right * std::log(spat_right));
  429 + }
  430 + }
  431 +
  432 + // Multiply-in normalisation
  433 + double spat = spat_raw * m_norm.value();
  434 +
  435 + // Make sure that spatial component does not become negative
  436 + if (spat < 0.0) {
  437 + spat = 0.0;
  438 + }
  439 +
  440 + // Divide spatial model value by event size
  441 + if (event.size() != 0.0) {
  442 + spat /= event.size();
  443 + }
  444 + else {
  445 + spat = 0.0;
  446 + }
  447 +
  448 + // Evaluate spectral and temporal components
  449 + double spec = (spectral() != NULL)
  450 + ? spectral()->eval(event.energy(), event.time(), gradients)
  451 + : 1.0;
  452 + double temp = (temporal() != NULL)
  453 + ? temporal()->eval(event.time(), gradients) : 1.0;
  454 +
  455 + // Compute value
  456 + double value = spat * spec * temp;
  457 +
  458 + // If gradients were requested then multiply factors to spectral and
  459 + // temporal gradients
  460 + if (gradients) {
  461 +
  462 + // Multiply factors to spatial gradient
  463 + if (m_norm.is_free()) {
  464 + double fact = spec * temp;
  465 + double g_norm = spat_raw * m_norm.scale();
  466 + m_norm.factor_gradient(g_norm * fact);
  467 + }
  468 +
  469 + // Multiply factors to spectral gradients
  470 + if (spectral() != NULL) {
  471 + double fact = spat * temp;
  472 + if (fact != 1.0) {
  473 + for (int i = 0; i < spectral()->size(); ++i)
  474 + (*spectral())[i].factor_gradient((*spectral())[i].factor_gradient() * fact);
  475 + }
  476 + }
  477 +
  478 + // Multiply factors to temporal gradients
  479 + if (temporal() != NULL) {
  480 + double fact = spat * spec;
  481 + if (fact != 1.0) {
  482 + for (int i = 0; i < temporal()->size(); ++i)
  483 + (*temporal())[i].factor_gradient((*temporal())[i].factor_gradient() * fact);
  484 + }
  485 + }
  486 +
  487 + } // endif: gradients were requested
  488 +
  489 + // Return value
  490 + return value;
  491 +}
  492 +
  493 +
  494 +/***********************************************************************//**
  495 + * @brief Return spatially integrated event rate in units of
  496 + * events MeV\f$^{-1}\f$ s\f$^{-1}\f$
  497 + *
  498 + * @param[in] obsEng Measured event energy.
  499 + * @param[in] obsTime Measured event time.
  500 + * @param[in] obs Observation.
  501 + * @return Spatially integrated event rate
  502 + * (events MeV\f$^{-1}\f$ s\f$^{-1}\f$)
  503 + *
  504 + * Spatially integrates the sky cube model for a given measured event energy
  505 + * and event time.
  506 + ***************************************************************************/
  507 +double GCTAModelSkyCube::npred(const GEnergy& obsEng,
  508 + const GTime& obsTime,
  509 + const GObservation& obs) const
  510 +{
  511 + // Initialise result
  512 + double npred = 0.0;
  513 + bool has_npred = false;
  514 +
  515 + // Build unique identifier
  516 + std::string id = obs.instrument() + ":" + obs.id();
  517 +
  518 + // Check if Npred value is already in cache
  519 + #if defined(G_USE_NPRED_CACHE)
  520 + if (!m_npred_names.empty()) {
  521 +
  522 + // Search for unique identifier, and if found, recover Npred value
  523 + // and break
  524 + for (int i = 0; i < m_npred_names.size(); ++i) {
  525 + if (m_npred_names[i] == id && m_npred_energies[i] == obsEng) {
  526 + npred = m_npred_values[i];
  527 + has_npred = true;
  528 + #if defined(G_DEBUG_NPRED)
  529 + std::cout << "GCTAModelSkyCube::npred:";
  530 + std::cout << " cache=" << i;
  531 + std::cout << " npred=" << npred << std::endl;
  532 + #endif
  533 + break;
  534 + }
  535 + }
  536 +
  537 + } // endif: there were values in the Npred cache
  538 + #endif
  539 +
  540 + // Continue only if no Npred cache value has been found
  541 + if (!has_npred) {
  542 +
  543 + // Evaluate only if model is valid
  544 + if (valid_model()) {
  545 +
  546 + // Set indices and weighting factors for energy interpolation
  547 + m_elogmeans.set_value(obsEng.log10TeV());
  548 + int inx_left = m_elogmeans.inx_left();
  549 + int inx_right = m_elogmeans.inx_right();
  550 + double wgt_left = m_elogmeans.wgt_left();
  551 + double wgt_right = m_elogmeans.wgt_right();
  552 +
  553 + // Loop over all map pixels
  554 + for (int i = 0; i < m_cube.npix(); ++i) {
  555 +
  556 + // Get bin value
  557 + double value = wgt_left * m_cube(i, inx_left) +
  558 + wgt_right * m_cube(i, inx_right);
  559 +
  560 + // Sum bin contents
  561 + npred += value * m_cube.solidangle(i);
  562 +
  563 + }
  564 +
  565 + // Store result in Npred cache
  566 + #if defined(G_USE_NPRED_CACHE)
  567 + m_npred_names.push_back(id);
  568 + m_npred_energies.push_back(obsEng);
  569 + m_npred_times.push_back(obsTime);
  570 + m_npred_values.push_back(npred);
  571 + #endif
  572 +
  573 + // Debug: Check for NaN
  574 + #if defined(G_NAN_CHECK)
  575 + if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
  576 + std::string origin = "GCTAModelSkyCube::npred";
  577 + std::string message = " NaN/Inf encountered (npred=" +
  578 + gammalib::str(npred) + ")";
  579 + gammalib::warning(origin, message);
  580 + }
  581 + #endif
  582 +
  583 + } // endif: model was valid
  584 +
  585 + } // endif: Npred computation required
  586 +
  587 + // Multiply in normalisation value
  588 + npred *= m_norm.value();
  589 +
  590 + // Multiply in spectral and temporal components
  591 + npred *= spectral()->eval(obsEng, obsTime);
  592 + npred *= temporal()->eval(obsTime);
  593 +
  594 + // Return Npred
  595 + return npred;
  596 +}
  597 +
  598 +
  599 +/***********************************************************************//**
  600 + * @brief Return simulated list of events
  601 + *
  602 + * @param[in] obs Observation.
  603 + * @param[in] ran Random number generator.
  604 + * @return Pointer to list of simulated events (needs to be de-allocated by
  605 + * client)
  606 + *
  607 + * @exception GException::feature_not_implemented
  608 + * Class does not support the generation of event lists.
  609 + *
  610 + * The simulation of an event list from a sky cube model is not implemented,
  611 + * hence the method will always throw an exception.
  612 + ***************************************************************************/
  613 +GCTAEventList* GCTAModelSkyCube::mc(const GObservation& obs, GRan& ran) const
  614 +{
  615 + // Feature not yet implemented
  616 + throw GException::feature_not_implemented(G_MC,
  617 + "MC computation not implemented for binned analysis.");
  618 +
  619 + // Return NULL pointer
  620 + return NULL;
  621 +
  622 +}
  623 +
  624 +
  625 +/***********************************************************************//**
  626 + * @brief Read CTA sky cube model from XML element
  627 + *
  628 + * @param[in] xml XML element.
  629 + *
  630 + * Set up CTA sky cube model from the information provided by an XML element.
  631 + * The XML element is expected to have the following structure
  632 + *
  633 + * <source name="CTA sky cube" type="CTASkyCube" instrument="CTA">
  634 + * <spatialModel type="ModelCube" file="model_cube.fits">
  635 + * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
  636 + * </spatialModel>
  637 + * <spectrum type="...">
  638 + * ...
  639 + * </spectrum>
  640 + * </source>
  641 + *
  642 + * Optionally, a temporal model may be provided using the following
  643 + * syntax
  644 + *
  645 + * <source name="CTA sky cube" type="CTASkyCube" instrument="CTA">
  646 + * <spatialModel type="ModelCube" file="model_cube.fits">
  647 + * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
  648 + * </spatialModel>
  649 + * <spectrum type="...">
  650 + * ...
  651 + * </spectrum>
  652 + * <temporal type="...">
  653 + * ...
  654 + * </temporal>
  655 + * </source>
  656 + *
  657 + * If no temporal component is found a constant model is assumed.
  658 + ***************************************************************************/
  659 +void GCTAModelSkyCube::read(const GXmlElement& xml)
  660 +{
  661 + // Clear model
  662 + clear();
  663 +
  664 + // Get pointers on spatial and spectral components
  665 + const GXmlElement* spat = xml.element("spatialModel", 0);
  666 + const GXmlElement* spec = xml.element("spectrum", 0);
  667 +
  668 + // Set spatial and spectral models
  669 + read_xml_spatial(*spat);
  670 + m_spectral = xml_spectral(*spec);
  671 +
  672 + // Handle optional temporal model
  673 + if (xml.elements("temporal")) {
  674 + const GXmlElement* temporal = xml.element("temporal", 0);
  675 + m_temporal = xml_temporal(*temporal);
  676 + }
  677 + else {
  678 + GModelTemporalConst constant;
  679 + m_temporal = constant.clone();
  680 + }
  681 +
  682 + // Read model attributes
  683 + read_attributes(xml);
  684 +
  685 + // Set parameter pointers
  686 + set_pointers();
  687 +
  688 + // Return
  689 + return;
  690 +}
  691 +
  692 +
  693 +/***********************************************************************//**
  694 + * @brief Write CTA cube background model into XML element
  695 + *
  696 + * @param[in] xml XML element.
  697 + *
  698 + * Write CTA cube background model information into an XML element.
  699 + * The XML element will have the following structure
  700 + *
  701 + * <source name="CTA sky cube" type="CTASkyCube" instrument="CTA">
  702 + * <spatialModel type="ModelCube" file="model_cube.fits">
  703 + * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
  704 + * </spatialModel>
  705 + * <spectrum type="...">
  706 + * ...
  707 + * </spectrum>
  708 + * </source>
  709 + *
  710 + * If the model contains a non-constant temporal model, the temporal
  711 + * component will also be written following the syntax
  712 + *
  713 + * <source name="CTA sky cube" type="CTASkyCube" instrument="CTA">
  714 + * <spatialModel type="ModelCube" file="model_cube.fits">
  715 + * <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
  716 + * </spatialModel>
  717 + * <spectrum type="...">
  718 + * ...
  719 + * </spectrum>
  720 + * <temporal type="...">
  721 + * ...
  722 + * </temporal>
  723 + * </source>
  724 + *
  725 + * If no temporal component is found a constant model is assumed.
  726 + ***************************************************************************/
  727 +void GCTAModelSkyCube::write(GXmlElement& xml) const
  728 +{
  729 + // Initialise pointer on source
  730 + GXmlElement* src = NULL;
  731 +
  732 + // Search corresponding source
  733 + int n = xml.elements("source");
  734 + for (int k = 0; k < n; ++k) {
  735 + GXmlElement* element = xml.element("source", k);
  736 + if (element->attribute("name") == name()) {
  737 + src = element;
  738 + break;
  739 + }
  740 + }
  741 +
  742 + // If we have a temporal model that is either not a constant, or a
  743 + // constant with a normalization value that differs from 1.0 then
  744 + // write the temporal component into the XML element. This logic
  745 + // assures compatibility with the Fermi/LAT format as this format
  746 + // does not handle temporal components.
  747 + bool write_temporal = ((m_temporal != NULL) &&
  748 + (m_temporal->type() != "Constant" ||
  749 + (*m_temporal)[0].value() != 1.0));
  750 +
  751 + // If no source with corresponding name was found then append one
  752 + if (src == NULL) {
  753 + src = xml.append("source");
  754 + src->append(GXmlElement("spatialModel"));
  755 + src->append(GXmlElement("spectrum"));
  756 + if (write_temporal) {
  757 + src->append(GXmlElement("temporal"));
  758 + }
  759 + }
  760 +
  761 + // Write spatial model
  762 + GXmlElement* spat = src->element("spatialModel", 0);
  763 + write_xml_spatial(*spat);
  764 +
  765 + // Write spectral model
  766 + if (spectral() != NULL) {
  767 + GXmlElement* spec = src->element("spectrum", 0);
  768 + spectral()->write(*spec);
  769 + }
  770 +
  771 + // Optionally write temporal model
  772 + if (write_temporal) {
  773 + if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
  774 + GXmlElement* temp = src->element("temporal", 0);
  775 + temporal()->write(*temp);
  776 + }
  777 + }
  778 +
  779 + // Write model attributes
  780 + write_attributes(*src);
  781 +
  782 + // Return
  783 + return;
  784 +}
  785 +
  786 +
  787 +/***********************************************************************//**
  788 + * @brief Print CTA cube background model information
  789 + *
  790 + * @param[in] chatter Chattiness.
  791 + * @return String containing CTA cube background model information.
  792 + ***************************************************************************/
  793 +std::string GCTAModelSkyCube::print(const GChatter& chatter) const
  794 +{
  795 + // Initialise result string
  796 + std::string result;
  797 +
  798 + // Continue only if chatter is not silent
  799 + if (chatter != SILENT) {
  800 +
  801 + // Append header
  802 + result.append("=== GCTAModelSkyCube ===");
  803 +
  804 + // Determine number of parameters per type
  805 + int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
  806 + int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
  807 +
  808 + // Append attributes
  809 + result.append("\n"+print_attributes());
  810 +
  811 + // Append model type
  812 + result.append("\n"+gammalib::parformat("Model type"));
  813 + if (n_spectral > 0) {
  814 + result.append("\""+spectral()->type()+"\"");
  815 + if (n_temporal > 0) {
  816 + result.append(" * ");
  817 + }
  818 + }
  819 + if (n_temporal > 0) {
  820 + result.append("\""+temporal()->type()+"\"");
  821 + }
  822 +
  823 + // Append parameters
  824 + result.append("\n"+gammalib::parformat("Number of parameters") +
  825 + gammalib::str(size()));
  826 + result.append("\n"+gammalib::parformat("Number of spectral par's") +
  827 + gammalib::str(n_spectral));
  828 + for (int i = 0; i < n_spectral; ++i) {
  829 + result.append("\n"+(*spectral())[i].print());
  830 + }
  831 + result.append("\n"+gammalib::parformat("Number of temporal par's") +
  832 + gammalib::str(n_temporal));
  833 + for (int i = 0; i < n_temporal; ++i) {
  834 + result.append("\n"+(*temporal())[i].print());
  835 + }
  836 +
  837 + } // endif: chatter was not silent
  838 +
  839 + // Return result
  840 + return result;
  841 +}
  842 +
  843 +
  844 +/*==========================================================================
  845 + = =
  846 + = Private methods =
  847 + = =
  848 + ==========================================================================*/
  849 +
  850 +/***********************************************************************//**
  851 + * @brief Initialise class members
  852 + ***************************************************************************/
  853 +void GCTAModelSkyCube::init_members(void)
  854 +{
  855 + // Initialise Value
  856 + m_norm.clear();
  857 + m_norm.name("Normalization");
  858 + m_norm.value(1.0);
  859 + m_norm.scale(1.0);
  860 + m_norm.range(0.001, 1000.0);
  861 + m_norm.gradient(0.0);
  862 + m_norm.fix();
  863 + m_norm.has_grad(true);
  864 +
  865 + // Set parameter pointer(s)
  866 + m_pars.clear();
  867 + m_pars.push_back(&m_norm);
  868 +
  869 + // Initialise members
  870 + m_filename.clear();
  871 + m_cube.clear();
  872 + m_ebounds.clear();
  873 + m_elogmeans.clear();
  874 + m_spectral = NULL;
  875 + m_temporal = NULL;
  876 +
  877 + // Initialise Npred cache
  878 + m_npred_names.clear();
  879 + m_npred_energies.clear();
  880 + m_npred_times.clear();
  881 + m_npred_values.clear();
  882 +
  883 + // Return
  884 + return;
  885 +}
  886 +
  887 +
  888 +/***********************************************************************//**
  889 + * @brief Copy class members
  890 + *
  891 + * @param[in] sky CTA sky cube model.
  892 + ***************************************************************************/
  893 +void GCTAModelSkyCube::copy_members(const GCTAModelSkyCube& sky)
  894 +{
  895 + // Copy members
  896 + m_filename = sky.m_filename;
  897 + m_cube = sky.m_cube;
  898 + m_ebounds = sky.m_ebounds;
  899 + m_elogmeans = sky.m_elogmeans;
  900 + m_norm = sky.m_norm;
  901 +
  902 + // Copy cache
  903 + m_npred_names = sky.m_npred_names;
  904 + m_npred_energies = sky.m_npred_energies;
  905 + m_npred_times = sky.m_npred_times;
  906 + m_npred_values = sky.m_npred_values;
  907 +
  908 + // Clone spectral and temporal model components
  909 + m_spectral = (sky.m_spectral != NULL) ? sky.m_spectral->clone() : NULL;
  910 + m_temporal = (sky.m_temporal != NULL) ? sky.m_temporal->clone() : NULL;
  911 +
  912 + // Set parameter pointers
  913 + set_pointers();
  914 +
  915 + // Return
  916 + return;
  917 +}
  918 +
  919 +
  920 +/***********************************************************************//**
  921 + * @brief Delete class members
  922 + ***************************************************************************/
  923 +void GCTAModelSkyCube::free_members(void)
  924 +{
  925 + // Free memory
  926 + if (m_spectral != NULL) delete m_spectral;
  927 + if (m_temporal != NULL) delete m_temporal;
  928 +
  929 + // Signal free pointers
  930 + m_spectral = NULL;
  931 + m_temporal = NULL;
  932 +
  933 + // Return
  934 + return;
  935 +}
  936 +
  937 +
  938 +/***********************************************************************//**
  939 + * @brief Set pointers
  940 + *
  941 + * Set pointers to all model parameters. The pointers are stored in a vector
  942 + * that is member of the GModelData base class.
  943 + ***************************************************************************/
  944 +void GCTAModelSkyCube::set_pointers(void)
  945 +{
  946 + // Clear parameters
  947 + m_pars.clear();
  948 +
  949 + // Push normalisation parameter on stack
  950 + m_pars.push_back(&m_norm);
  951 +
  952 + // Determine the number of parameters
  953 + int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
  954 + int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
  955 + int n_pars = n_spectral + n_temporal;
  956 +
  957 + // Continue only if there are parameters
  958 + if (n_pars > 0) {
  959 +
  960 + // Gather spectral parameters
  961 + for (int i = 0; i < n_spectral; ++i) {
  962 + m_pars.push_back(&((*spectral())[i]));
  963 + }
  964 +
  965 + // Gather temporal parameters
  966 + for (int i = 0; i < n_temporal; ++i) {
  967 + m_pars.push_back(&((*temporal())[i]));
  968 + }
  969 +
  970 + }
  971 +
  972 + // Return
  973 + return;
  974 +}
  975 +
  976 +
  977 +/***********************************************************************//**
  978 + * @brief Verifies if model has all components
  979 + *
  980 + * Returns 'true' if models has a spectral and a temporal component.
  981 + * Otherwise returns 'false'.
  982 + ***************************************************************************/
  983 +bool GCTAModelSkyCube::valid_model(void) const
  984 +{
  985 + // Set result
  986 + bool result = ((spectral() != NULL) && (temporal() != NULL));
  987 +
  988 + // Return result
  989 + return result;
  990 +}
  991 +
  992 +
  993 +/***********************************************************************//**
  994 + * @brief Read spatial model from XML element
  995 + *
  996 + * @param[in] xml XML element.
  997 + *
  998 + * @exception GException::invalid_value
  999 + * Invalid XML format encountered.
  1000 + *
  1001 + * Read sky cube information from the spatial XML element. The expected
  1002 + * format of the XML element is
  1003 + *
  1004 + * <spatialModel type="ModelCube" file="model_cube.fits">
  1005 + * <parameter name="Normalization" ../>
  1006 + * </spatialModel>
  1007 + ***************************************************************************/
  1008 +void GCTAModelSkyCube::read_xml_spatial(const GXmlElement& xml)
  1009 +{
  1010 + // Verify model type
  1011 + if (xml.attribute("type") != "ModelCube") {
  1012 + std::string msg = "Spatial model type \""+xml.attribute("type")+
  1013 + "\" is not of type \"ModelCube\". Please verify "
  1014 + "the XML format.";
  1015 + throw GException::invalid_value(G_READ_XML_SPATIAL, msg);
  1016 + }
  1017 +
  1018 + // Get filename
  1019 + GFilename filename = gammalib::xml_file_expand(xml, xml.attribute("file"));
  1020 +
  1021 + // Get XML parameter
  1022 + const GXmlElement* norm = gammalib::xml_get_par(G_READ_XML_SPATIAL, xml,
  1023 + m_norm.name());
  1024 +
  1025 + // Read parameter
  1026 + m_norm.read(*norm);
  1027 +
  1028 + // Load sky cube
  1029 + load(filename);
  1030 +
  1031 + // Return
  1032 + return;
  1033 +}
  1034 +
  1035 +
  1036 +/***********************************************************************//**
  1037 + * @brief Write spatial model into XML element
  1038 + *
  1039 + * @param[in] xml XML element.
  1040 + *
  1041 + * @exception GException::invalid_value
  1042 + * Invalid XML format encountered.
  1043 + *
  1044 + * Write sky cube information into the spatial XML element. The format of
  1045 + * the XML element is
  1046 + *
  1047 + * <spatialModel type="ModelCube" file="model_cube.fits">
  1048 + * <parameter name="Normalization" ../>
  1049 + * </spatialModel>
  1050 + ***************************************************************************/
  1051 +void GCTAModelSkyCube::write_xml_spatial(GXmlElement& xml) const
  1052 +{
  1053 + // Set model type
  1054 + if (xml.attribute("type") == "") {
  1055 + xml.attribute("type", "ModelCube");
  1056 + }
  1057 +
  1058 + // Verify model type
  1059 + if (xml.attribute("type") != "ModelCube") {
  1060 + std::string msg = "Spatial model type \""+xml.attribute("type")+
  1061 + "\" is not of type \"ModelCube\". Please verify "
  1062 + "the XML format.";
  1063 + throw GException::invalid_value(G_WRITE_XML_SPATIAL, msg);
  1064 + }
  1065 +
  1066 + // Set filename
  1067 + xml.attribute("file", gammalib::xml_file_reduce(xml, m_filename));
  1068 +
  1069 + // Get XML parameter
  1070 + GXmlElement* norm = gammalib::xml_need_par(G_WRITE_XML_SPATIAL, xml,
  1071 + m_norm.name());
  1072 +
  1073 + // Write parameter
  1074 + m_norm.write(*norm);
  1075 +
  1076 + // Return
  1077 + return;
  1078 +}
  1079 +
  1080 +
  1081 +/***********************************************************************//**
  1082 + * @brief Return pointer to spectral model from XML element
  1083 + *
  1084 + * @param[in] spectral XML element.
  1085 + * @return Pointer to spectral model.
  1086 + *
  1087 + * Returns pointer to spectral model that is defined in an XML element.
  1088 + ***************************************************************************/
  1089 +GModelSpectral* GCTAModelSkyCube::xml_spectral(const GXmlElement& spectral) const
  1090 +{
  1091 + // Get spectral model
  1092 + GModelSpectralRegistry registry;
  1093 + GModelSpectral* ptr = registry.alloc(spectral);
  1094 +
  1095 + // Return pointer
  1096 + return ptr;
  1097 +}
  1098 +
  1099 +
  1100 +/***********************************************************************//**
  1101 + * @brief Return pointer to temporal model from XML element
  1102 + *
  1103 + * @param[in] temporal XML element.
  1104 + * @return Pointer to temporal model.
  1105 + *
  1106 + * Returns pointer to temporal model that is defined in an XML element.
  1107 + ***************************************************************************/
  1108 +GModelTemporal* GCTAModelSkyCube::xml_temporal(const GXmlElement& temporal) const
  1109 +{
  1110 + // Get temporal model
  1111 + GModelTemporalRegistry registry;
  1112 + GModelTemporal* ptr = registry.alloc(temporal);
  1113 +
  1114 + // Return pointer
  1115 + return ptr;
  1116 +}
... ...
inst/cta/test/data/crab_modcube.fits.gz 0 → 100644
View file @4f9fa41
No preview for this file type
inst/cta/test/data/cta_model_sky_cube.xml 0 → 100644
  View file @4f9fa41
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="CTA sky cube" type="CTASkyCube" instrument="CTA">
  4 + <spatialModel type="ModelCube" file="crab_modcube.fits">
  5 + <parameter name="Normalization" scale="1" value="1" min="0.1" max="10" free="0"/>
  6 + </spatialModel>
  7 + <spectrum type="PowerLaw">
  8 + <parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
  9 + <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
  10 + <parameter name="PivotEnergy" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
  11 + </spectrum>
  12 + </source>
  13 +</source_library>
... ...
inst/cta/test/test_CTA.cpp
  View file @4f9fa41
... ... @@ -54,10 +54,12 @@ const std::string cta_irf = &quot;cta_dummy_irf&quot;;
54 54 const std::string cta_events = datadir+"/crab_events.fits";
55 55 const std::string cta_events_gti = datadir+"/crab_events_gti.fits[EVENTS2]";
56 56 const std::string cta_cntmap = datadir+"/crab_cntmap.fits";
  57 +const std::string cta_modcube = datadir+"/crab_modcube.fits";
57 58 const std::string cta_bin_xml = datadir+"/obs_binned.xml";
58 59 const std::string cta_unbin_xml = datadir+"/obs_unbinned.xml";
59 60 const std::string cta_model_xml = datadir+"/crab.xml";
60 61 const std::string cta_rsp_xml = datadir+"/rsp_models.xml";
  62 +const std::string cta_sky_cube_xml = datadir+"/cta_model_sky_cube.xml";
61 63 const std::string cta_bgd_rad_gauss_xml = datadir+"/cta_model_bgd_rad_gauss.xml";
62 64 const std::string cta_bgd_gauss_e_xml = datadir+"/cta_model_bgd_gauss_e.xml";
63 65 const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml";
... ... @@ -171,6 +173,8 @@ void TestGCTAModel::set(void)
171 173 name("Test CTA models");
172 174  
173 175 // Append tests to test suite
  176 + append(static_cast<pfunction>(&TestGCTAModel::test_model_sky_cube),
  177 + "Test CTA sky cube model");
174 178 append(static_cast<pfunction>(&TestGCTAModel::test_model_bgd),
175 179 "Test CTA background model");
176 180 append(static_cast<pfunction>(&TestGCTAModel::test_model_cube_bgd),
... ... @@ -1748,6 +1752,67 @@ void TestGCTAResponse::test_edispcube_integration(const GCTACubeEdisp&amp; edisp,
1748 1752  
1749 1753  
1750 1754 /***********************************************************************//**
  1755 + * @brief Test CTA sky cube model
  1756 + ***************************************************************************/
  1757 +void TestGCTAModel::test_model_sky_cube(void)
  1758 +{
  1759 + // Test void constuctor
  1760 + GCTAModelSkyCube model1;
  1761 + test_value(model1.type(), "CTASkyCube", "Check void model type");
  1762 + test_value(model1.size(), 1, "Check void model size");
  1763 + test_value(model1["Normalization"].value(), 1.0, "Check void model \"Normalization\" value");
  1764 + test_assert(model1.is_constant(), "Check if void model is constant.");
  1765 +
  1766 + // Test XML constructor
  1767 + GXml xml(cta_sky_cube_xml);
  1768 + const GXmlElement& lib = *xml.element("source_library", 0);
  1769 + const GXmlElement& src = *lib.element("source", 0);
  1770 + GCTAModelSkyCube model2(src);
  1771 + test_value(model2.size(), 5, "Check model2 size");
  1772 + test_value(model2["Normalization"].value(), 1.0, "Check model2 \"Normalization\" value");
  1773 + test_value(model2["Prefactor"].value(), 5.7e-16, "Check model2 \"Prefactor\" value");
  1774 + test_value(model2["Index"].value(), -2.48, "Check model2 \"Index\" value");
  1775 + test_value(model2["PivotEnergy"].value(), 0.3e6, "Check model2 \"PivotEnergy\" value");
  1776 + test_value(model2.filename(), cta_modcube, "Check model2 filename");
  1777 + test_assert(model2.is_constant(), "Check if model2 is constant.");
  1778 +
  1779 + // Test filename and spectral constuctor
  1780 + GModelSpectralPlaw plaw(1.0, 0.0, GEnergy(1.0, "TeV"));
  1781 + GCTAModelSkyCube model3(cta_modcube, plaw);
  1782 + test_value(model3.size(), 5, "Check model3 size");
  1783 + test_value(model3["Normalization"].value(), 1.0, "Check model3 \"Normalization\" value");
  1784 + test_value(model3["Prefactor"].value(), 1.0, "Check model3 \"Prefactor\" value");
  1785 + test_value(model3["Index"].value(), 0.0, "Check model3 \"Index\" value");
  1786 + test_value(model3["PivotEnergy"].value(), 1.0e6, "Check model3 \"PivotEnergy\" value");
  1787 + test_value(model3.filename(), cta_modcube, "Check model3 filename");
  1788 + test_assert(model3.is_constant(), "Check if model3 is constant.");
  1789 +
  1790 + // Test loading of sky cube in model container
  1791 + GModels models(cta_sky_cube_xml);
  1792 + GModel* model = models["CTA sky cube"];
  1793 + test_value(model->size(), 5, "Check model size");
  1794 + test_value((*model)["Normalization"].value(), 1.0, "Check model \"Normalization\" value");
  1795 + test_value((*model)["Prefactor"].value(), 5.7e-16, "Check model \"Prefactor\" value");
  1796 + test_value((*model)["Index"].value(), -2.48, "Check model \"Index\" value");
  1797 + test_value((*model)["PivotEnergy"].value(), 0.3e6, "Check model \"PivotEnergy\" value");
  1798 + test_assert(model->is_constant(), "Check if model is constant.");
  1799 +
  1800 + // Test XML saving and reloading of cube background
  1801 + models.save("test.xml");
  1802 + models.load("test.xml");
  1803 + model = models["CTA sky cube"];
  1804 + test_value(model->size(), 5, "Check model size");
  1805 + test_value((*model)["Normalization"].value(), 1.0, "Check model \"Normalization\" value");
  1806 + test_value((*model)["Prefactor"].value(), 5.7e-16, "Check model \"Prefactor\" value");
  1807 + test_value((*model)["Index"].value(), -2.48, "Check model \"Index\" value");
  1808 + test_value((*model)["PivotEnergy"].value(), 0.3e6, "Check model \"PivotEnergy\" value");
  1809 + test_assert(model->is_constant(), "Check if model is constant.");
  1810 +
  1811 + // Return
  1812 + return;
  1813 +
  1814 +}
  1815 +/***********************************************************************//**
1751 1816 * @brief Test CTA background
1752 1817 ***************************************************************************/
1753 1818 void TestGCTAModel::test_model_bgd(void)
... ...
inst/cta/test/test_CTA.hpp
  View file @4f9fa41
1 1 /***************************************************************************
2 2 * test_CTA.hpp - Test CTA classes *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -121,6 +121,7 @@ public:
121 121 virtual void set(void);
122 122 virtual TestGCTAModel* clone(void) const { return new TestGCTAModel(*this); }
123 123 virtual std::string classname(void) const { return "TestGCTAModel"; }
  124 + void test_model_sky_cube(void);
124 125 void test_model_bgd(void);
125 126 void test_model_cube_bgd(void);
126 127 void test_model_irf_bgd(void);
... ...
inst/cta/test/test_CTA.py
  View file @4f9fa41
1 1 # ==========================================================================
2 2 # This module performs unit tests for the GammaLib CTA module.
3 3 #
4   -# Copyright (C) 2012-2018 Juergen Knoedlseder
  4 +# Copyright (C) 2012-2020 Juergen Knoedlseder
5 5 #
6 6 # This program is free software: you can redistribute it and/or modify
7 7 # it under the terms of the GNU General Public License as published by
... ... @@ -405,6 +405,7 @@ class Test(gammalib.GPythonTestSuite):
405 405 test_support.pickeling(self, gammalib.GCTAEventCube())
406 406 test_support.pickeling(self, gammalib.GCTAEventList())
407 407 test_support.pickeling(self, gammalib.GCTAInstDir())
  408 + test_support.pickeling(self, gammalib.GCTAModelSkyCube())
408 409 test_support.pickeling(self, gammalib.GCTAModelBackground())
409 410 test_support.pickeling(self, gammalib.GCTAModelAeffBackground())
410 411 test_support.pickeling(self, gammalib.GCTAModelCubeBackground())
... ... @@ -491,6 +492,7 @@ class Test(gammalib.GPythonTestSuite):
491 492 hessirf.load_psf(hesspsf)
492 493 hessirf.load_edisp(hessedisp)
493 494 hessirf.load_background(hessbkg)
  495 + skycube = self._data + '/crab_modcube.fits.gz'
494 496  
495 497 # Perform pickeling tests of filled classes
496 498 test_support.pickeling(self, gammalib.GCTAAeff2D(irfname))
... ... @@ -512,6 +514,7 @@ class Test(gammalib.GPythonTestSuite):
512 514 test_support.pickeling(self, gammalib.GCTAEventCube(cube))
513 515 test_support.pickeling(self, gammalib.GCTAEventList(list))
514 516 test_support.pickeling(self, gammalib.GCTAInstDir(instdir))
  517 + test_support.pickeling(self, gammalib.GCTAModelSkyCube(skycube,plaw))
515 518 test_support.pickeling(self, gammalib.GCTAModelBackground(radgauss,plaw))
516 519 test_support.pickeling(self, gammalib.GCTAModelAeffBackground(plaw))
517 520 test_support.pickeling(self, gammalib.GCTAModelCubeBackground(plaw))
... ...