Sergio Hernández Cadena / gammalib

Commit caa5b1ed7d203160c4e83fe5ebb47e68bd67ba53

Authored by Michael Mayer
1 parent 86b32056

Added GCTAModelAeffBackground class and introduced GCTAAeff::max function

Inline Side-by-side
Showing 20 changed files with 1848 additions and 0 deletions
inst/cta/Makefile.am
  View file @caa5b1e
... ... @@ -49,6 +49,7 @@ sources = src/GCTASupport.cpp \
49 49 src/GCTAPointing.cpp \
50 50 src/GCTAModelCubeBackground.cpp \
51 51 src/GCTAModelIrfBackground.cpp \
  52 + src/GCTAModelAeffBackground.cpp \
52 53 src/GCTAModelRadialRegistry.cpp \
53 54 src/GCTAModelRadial.cpp \
54 55 src/GCTAModelRadialGauss.cpp \
... ... @@ -97,6 +98,7 @@ pkginclude_HEADERS = include/GCTATypemaps.hpp \
97 98 include/GCTACubeSourceDiffuse.hpp \
98 99 include/GCTAModelCubeBackground.hpp \
99 100 include/GCTAModelIrfBackground.hpp \
  101 + include/GCTAModelAeffBackground.hpp \
100 102 include/GCTAModelRadialRegistry.hpp \
101 103 include/GCTAModelRadial.hpp \
102 104 include/GCTAModelRadialGauss.hpp \
... ...
inst/cta/include/GCTAAeff.hpp
  View file @caa5b1e
... ... @@ -73,6 +73,10 @@ public:
73 73 virtual void load(const std::string& filename) = 0;
74 74 virtual std::string filename(void) const = 0;
75 75 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
  76 + virtual double max(const double& logE,
  77 + const double& zenith,
  78 + const double& azimuth,
  79 + const bool& etrue = true) const = 0;
76 80  
77 81 protected:
78 82 // Methods
... ...
inst/cta/include/GCTAAeff2D.hpp
  View file @caa5b1e
... ... @@ -67,6 +67,10 @@ public:
67 67 void load(const std::string& filename);
68 68 std::string filename(void) const;
69 69 std::string print(const GChatter& chatter = NORMAL) const;
  70 + double max(const double& logE,
  71 + const double& zenith,
  72 + const double& azimuth,
  73 + const bool& etrue = true) const;
70 74  
71 75 // Methods
72 76 const GCTAResponseTable& table(void) const;
... ...
inst/cta/include/GCTAAeffArf.hpp
  View file @caa5b1e
... ... @@ -72,6 +72,10 @@ public:
72 72 void load(const std::string& filename);
73 73 std::string filename(void) const;
74 74 std::string print(const GChatter& chatter = NORMAL) const;
  75 + double max(const double& logE,
  76 + const double& zenith,
  77 + const double& azimuth,
  78 + const bool& etrue = true) const;
75 79  
76 80 // Methods
77 81 int size(void) const;
... ...
inst/cta/include/GCTAAeffPerfTable.hpp
  View file @caa5b1e
... ... @@ -69,6 +69,10 @@ public:
69 69 void load(const std::string& filename);
70 70 std::string filename(void) const;
71 71 std::string print(const GChatter& chatter = NORMAL) const;
  72 + double max(const double& logE,
  73 + const double& zenith,
  74 + const double& azimuth,
  75 + const bool& etrue = true) const;
72 76  
73 77 // Methods
74 78 int size(void) const;
... ...
inst/cta/include/GCTALib.hpp
  View file @caa5b1e
... ... @@ -73,6 +73,7 @@
73 73 #include "GCTACubeSourceDiffuse.hpp"
74 74 #include "GCTAModelCubeBackground.hpp"
75 75 #include "GCTAModelIrfBackground.hpp"
  76 +#include "GCTAModelAeffBackground.hpp"
76 77 #include "GCTAModelRadial.hpp"
77 78 #include "GCTAModelRadialRegistry.hpp"
78 79 #include "GCTAModelRadialGauss.hpp"
... ...
inst/cta/include/GCTAModelAeffBackground.hpp 0 → 100644
  View file @caa5b1e
  1 +/***************************************************************************
  2 + * GCTAModelAeffBackground.hpp - CTA Aeff background model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 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 GCTAModelAeffBackground.hpp
  23 + * @brief CTA Aeff background model class definition
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +#ifndef GCTAMODELAEFFBACKGROUND_HPP
  28 +#define GCTAMODELAEFFBACKGROUND_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GModelData.hpp"
  33 +#include "GModelSpectral.hpp"
  34 +#include "GModelTemporal.hpp"
  35 +#include "GFunction.hpp"
  36 +#include "GCTAEventList.hpp"
  37 +#include "GCTAAeff.hpp"
  38 +
  39 +/* __ Forward declarations _______________________________________________ */
  40 +
  41 +
  42 +/***********************************************************************//**
  43 + * @class GCTAModelIrfBackground
  44 + *
  45 + * @brief CTA IRF background model class
  46 + ***************************************************************************/
  47 +class GCTAModelAeffBackground : public GModelData {
  48 +
  49 +public:
  50 + // Constructors and destructors
  51 + GCTAModelAeffBackground(void);
  52 + explicit GCTAModelAeffBackground(const GXmlElement& xml);
  53 + explicit GCTAModelAeffBackground(const GModelSpectral& spectral);
  54 + GCTAModelAeffBackground(const GCTAModelAeffBackground& bgd);
  55 + virtual ~GCTAModelAeffBackground(void);
  56 +
  57 + // Operators
  58 + GCTAModelAeffBackground& operator=(const GCTAModelAeffBackground& bgd);
  59 +
  60 + // Implemented pure virtual methods
  61 + virtual void clear(void);
  62 + virtual GCTAModelAeffBackground* clone(void) const;
  63 + virtual std::string classname(void) const;
  64 + virtual std::string type(void) const;
  65 + virtual bool is_constant(void) const;
  66 + virtual double eval(const GEvent& event,
  67 + const GObservation& obs) const;
  68 + virtual double eval_gradients(const GEvent& event,
  69 + const GObservation& obs) const;
  70 + virtual double npred(const GEnergy& obsEng,
  71 + const GTime& obsTime,
  72 + const GObservation& obs) const;
  73 + virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
  74 + virtual void read(const GXmlElement& xml);
  75 + virtual void write(GXmlElement& xml) const;
  76 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  77 +
  78 + // Other methods
  79 + GModelSpectral* spectral(void) const;
  80 + GModelTemporal* temporal(void) const;
  81 +
  82 +private:
  83 + // Methods
  84 + void init_members(void);
  85 + void copy_members(const GCTAModelAeffBackground& bgd);
  86 + void free_members(void);
  87 + void set_pointers(void);
  88 + double aeff_integral(const GObservation& obs, const double& logE) const;
  89 + bool valid_model(void) const;
  90 + GModelSpectral* xml_spectral(const GXmlElement& spectral) const;
  91 + GModelTemporal* xml_temporal(const GXmlElement& temporal) const;
  92 +
  93 + // ROI integration kernel over theta
  94 + class npred_roi_kern_theta : public GFunction {
  95 + public:
  96 + npred_roi_kern_theta(const GCTAAeff* aeff,
  97 + const double& logE,
  98 + const int& iter) :
  99 + m_aeff(aeff),
  100 + m_logE(logE),
  101 + m_iter(iter) { }
  102 + double eval(const double& theta);
  103 + protected:
  104 + const GCTAAeff* m_aeff; //!< Pointer to background
  105 + const double& m_logE; //!< Log10 of energy
  106 + const int& m_iter; //!< Romberg iterations
  107 + };
  108 +
  109 + // ROI integration kernel over phi
  110 + class npred_roi_kern_phi : public GFunction {
  111 + public:
  112 + npred_roi_kern_phi(const GCTAAeff* aeff,
  113 + const double& logE,
  114 + const double& theta) :
  115 + m_aeff(aeff),
  116 + m_logE(logE),
  117 + m_theta(theta) { }
  118 + double eval(const double& phi);
  119 + protected:
  120 + const GCTAAeff* m_aeff; //!< Pointer to background
  121 + const double& m_logE; //!< Log10 of energy
  122 + const double& m_theta; //!< Offset angle (radians)
  123 + };
  124 +
  125 + // Members
  126 + GModelSpectral* m_spectral; //!< Spectral model
  127 + GModelTemporal* m_temporal; //!< Temporal model
  128 +
  129 + // Npred cache
  130 + mutable std::vector<std::string> m_npred_names; //!< Model names
  131 + mutable std::vector<GEnergy> m_npred_energies; //!< Model energy
  132 + mutable std::vector<GTime> m_npred_times; //!< Model time
  133 + mutable std::vector<double> m_npred_values; //!< Model values
  134 +};
  135 +
  136 +
  137 +/***********************************************************************//**
  138 + * @brief Return class name
  139 + *
  140 + * @return String containing the class name ("GCTAModelAeffBackground").
  141 + ***************************************************************************/
  142 +inline
  143 +std::string GCTAModelAeffBackground::classname(void) const
  144 +{
  145 + return ("GCTAModelAeffBackground");
  146 +}
  147 +
  148 +
  149 +/***********************************************************************//**
  150 + * @brief Return data model type
  151 + *
  152 + * @return Data model type.
  153 + *
  154 + * Returns the type of the data model.
  155 + ***************************************************************************/
  156 +inline
  157 +std::string GCTAModelAeffBackground::type(void) const
  158 +{
  159 + return ("CTAAeffBackground");
  160 +}
  161 +
  162 +
  163 +/***********************************************************************//**
  164 + * @brief Signals if sky model is temporally constant
  165 + *
  166 + * @return True if sky model is temporally constant, false otherwise.
  167 + *
  168 + * Signals if the sky model is temporally constant. A temporally constant
  169 + * model is a model that has a temporal component of type "Constant".
  170 + ***************************************************************************/
  171 +inline
  172 +bool GCTAModelAeffBackground::is_constant(void) const
  173 +{
  174 + return (m_temporal != NULL && m_temporal->type() == "Constant");
  175 +}
  176 +
  177 +
  178 +/***********************************************************************//**
  179 + * @brief Return spectral model component
  180 + *
  181 + * @return Pointer to spectral model component.
  182 + *
  183 + * Returns a pointer to the spectral model component of the model. The
  184 + * pointer is of type GModelSpectral. Note that a NULL pointer may be
  185 + * returned if the sky model has no spectral model component.
  186 + ***************************************************************************/
  187 +inline
  188 +GModelSpectral* GCTAModelAeffBackground::spectral(void) const
  189 +{
  190 + return (m_spectral);
  191 +}
  192 +
  193 +
  194 +/***********************************************************************//**
  195 + * @brief Return temporal model component
  196 + *
  197 + * @return Pointer to temporal model component.
  198 + *
  199 + * Returns a pointer to the temporal model component of the model. The
  200 + * pointer is of type GModelTemporal. Note that a NULL pointer may be
  201 + * returned if the sky model has no temporal model component.
  202 + ***************************************************************************/
  203 +inline
  204 +GModelTemporal* GCTAModelAeffBackground::temporal(void) const
  205 +{
  206 + return (m_temporal);
  207 +}
  208 +
  209 +#endif /* GCTAMODELAEFFBACKGROUND_HPP */
... ...
inst/cta/pyext/GCTAAeff.i
  View file @caa5b1e
... ... @@ -57,6 +57,10 @@ public:
57 57 virtual std::string classname(void) const = 0;
58 58 virtual void load(const std::string& filename) = 0;
59 59 virtual std::string filename(void) const = 0;
  60 + virtual double max(const double& logE,
  61 + const double& zenith,
  62 + const double& azimuth,
  63 + const bool& etrue = true) const = 0;
60 64 };
61 65  
62 66  
... ...
inst/cta/pyext/GCTAAeff2D.i
  View file @caa5b1e
... ... @@ -58,6 +58,10 @@ public:
58 58 std::string classname(void) const;
59 59 void load(const std::string& filename);
60 60 std::string filename(void) const;
  61 + double max(const double& logE,
  62 + const double& zenith,
  63 + const double& azimuth,
  64 + const bool& etrue = true) const;
61 65  
62 66 // Methods
63 67 const GCTAResponseTable& table(void) const;
... ...
inst/cta/pyext/GCTAAeffArf.i
  View file @caa5b1e
... ... @@ -58,6 +58,10 @@ public:
58 58 std::string classname(void) const;
59 59 void load(const std::string& filename);
60 60 std::string filename(void) const;
  61 + double max(const double& logE,
  62 + const double& zenith,
  63 + const double& azimuth,
  64 + const bool& etrue = true) const;
61 65  
62 66 // Methods
63 67 int size(void) const;
... ...
inst/cta/pyext/GCTAAeffPerfTable.i
  View file @caa5b1e
... ... @@ -58,6 +58,10 @@ public:
58 58 std::string classname(void) const;
59 59 void load(const std::string& filename);
60 60 std::string filename(void) const;
  61 + double max(const double& logE,
  62 + const double& zenith,
  63 + const double& azimuth,
  64 + const bool& etrue = true) const;
61 65  
62 66 // Methods
63 67 int size(void) const;
... ...
inst/cta/pyext/GCTAModelAeffBackground.i 0 → 100644
  View file @caa5b1e
  1 +/***************************************************************************
  2 + * GCTAModelAeffBackground.i - CTA IRF background model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014-2015 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 GCTAModelAeffBackground.i
  23 + * @brief CTA IRF background model class definition
  24 + * @author Michael Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GCTAModelAeffBackground.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GCTAModelAeffBackground
  34 + *
  35 + * @brief CTA IRF background model class
  36 + ***************************************************************************/
  37 +class GCTAModelAeffBackground : public GModelData {
  38 +public:
  39 + // Constructors and destructors
  40 + GCTAModelAeffBackground(void);
  41 + explicit GCTAModelAeffBackground(const GXmlElement& xml);
  42 + explicit GCTAModelAeffBackground(const GModelSpectral& spectral);
  43 + GCTAModelAeffBackground(const GCTAModelAeffBackground& bgd);
  44 + virtual ~GCTAModelAeffBackground(void);
  45 +
  46 + // Implemented pure virtual methods
  47 + virtual void clear(void);
  48 + virtual GCTAModelAeffBackground* clone(void) const;
  49 + virtual std::string classname(void) const;
  50 + virtual std::string type(void) const;
  51 + virtual bool is_constant(void) const;
  52 + virtual double eval(const GEvent& event,
  53 + const GObservation& obs) const;
  54 + virtual double eval_gradients(const GEvent& event,
  55 + const GObservation& obs) const;
  56 + virtual double npred(const GEnergy& obsEng,
  57 + const GTime& obsTime,
  58 + const GObservation& obs) const;
  59 + virtual GCTAEventList* mc(const GObservation& obs, GRan& ran) const;
  60 + virtual void read(const GXmlElement& xml);
  61 + virtual void write(GXmlElement& xml) const;
  62 +
  63 + // Other methods
  64 + GModelSpectral* spectral(void) const;
  65 + GModelTemporal* temporal(void) const;
  66 +};
  67 +
  68 +/***********************************************************************//**
  69 +* @brief GCTAModelAeffBackground class extension
  70 +***************************************************************************/
  71 +%extend GCTAModelAeffBackground {
  72 + GCTAModelAeffBackground copy() {
  73 + return (*self);
  74 + }
  75 +};
... ...
inst/cta/pyext/cta.i
  View file @caa5b1e
... ... @@ -195,6 +195,7 @@
195 195 %include "GCTACubePsf.i"
196 196 %include "GCTAModelCubeBackground.i"
197 197 %include "GCTAModelIrfBackground.i"
  198 +%include "GCTAModelAeffBackground.i"
198 199 %include "GCTAModelRadial.i"
199 200 %include "GCTAModelRadialRegistry.i"
200 201 %include "GCTAModelRadialGauss.i"
... ...
inst/cta/src/GCTAAeff2D.cpp
  View file @caa5b1e
... ... @@ -408,6 +408,51 @@ std::string GCTAAeff2D::print(const GChatter&amp; chatter) const
408 408 }
409 409  
410 410  
  411 +/***********************************************************************//**
  412 + * @brief Return maximum effective area at a given energy
  413 + *
  414 + * @param[in] logE Log10 of the true photon energy (TeV).
  415 + * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
  416 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
  417 + * @param[in] etrue Use true energy (true/false). Not used.
  418 + * @return String containing effective area information.
  419 + ***************************************************************************/
  420 +double GCTAAeff2D::max(const double& logE,
  421 + const double& zenith,
  422 + const double& azimuth,
  423 + const bool& etrue) const
  424 +{
  425 + // Set parameter index
  426 + int index = (etrue) ? 0 : 1;
  427 +
  428 + int n_theta = m_aeff.axis(1);
  429 +
  430 + // Initialise maximum effective area
  431 + double max_aeff = 0.0;
  432 +
  433 + // Compute stepsize
  434 + double theta_step = m_aeff.axis_hi(1, n_theta - 1) / (double(n_theta) - 1.0);
  435 +
  436 + // Loop over theta values
  437 + for(int i = 0; i < m_aeff.axis(1); i++) {
  438 +
  439 + // Compute theta
  440 + double theta = i * theta_step;
  441 +
  442 + // Get effective area value in cm2
  443 + double aeff = m_aeff(index, logE, theta);
  444 +
  445 + // Adjust maximum effective area if necessary
  446 + if (aeff > max_aeff) {
  447 + max_aeff = aeff;
  448 + }
  449 + } // endfor: loop over theta values
  450 +
  451 + // Return effective area value
  452 + return max_aeff;
  453 +
  454 +}
  455 +
411 456 /*==========================================================================
412 457 = =
413 458 = Private methods =
... ...
inst/cta/src/GCTAAeffArf.cpp
  View file @caa5b1e
... ... @@ -477,6 +477,34 @@ std::string GCTAAeffArf::print(const GChatter&amp; chatter) const
477 477 }
478 478  
479 479  
  480 +/***********************************************************************//**
  481 + * @brief Return maximum effective area at a given energy
  482 + *
  483 + * @param[in] logE Log10 of the true photon energy (TeV).
  484 + * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
  485 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
  486 + * @param[in] etrue Use true energy (true/false). Not used.
  487 + * @return String containing effective area information.
  488 + ***************************************************************************/
  489 +double GCTAAeffArf::max(const double& logE,
  490 + const double& zenith,
  491 + const double& azimuth,
  492 + const bool& etrue) const
  493 +{
  494 +
  495 + // Get effective area value in cm2
  496 + double aeff_max = m_logE.interpolate(logE, m_aeff);
  497 +
  498 + // Make sure that effective area is not negative
  499 + if (aeff_max < 0.0) {
  500 + aeff_max = 0.0;
  501 + }
  502 +
  503 + // Return result
  504 + return aeff_max;
  505 +
  506 +}
  507 +
480 508 /*==========================================================================
481 509 = =
482 510 = Private methods =
... ...
inst/cta/src/GCTAAeffPerfTable.cpp
  View file @caa5b1e
... ... @@ -351,6 +351,34 @@ std::string GCTAAeffPerfTable::print(const GChatter&amp; chatter) const
351 351 }
352 352  
353 353  
  354 +/***********************************************************************//**
  355 + * @brief Return maximum effective area at a given energy
  356 + *
  357 + * @param[in] logE Log10 of the true photon energy (TeV).
  358 + * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
  359 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
  360 + * @param[in] etrue Use true energy (true/false). Not used.
  361 + * @return String containing effective area information.
  362 + ***************************************************************************/
  363 +double GCTAAeffPerfTable::max(const double& logE,
  364 + const double& zenith,
  365 + const double& azimuth,
  366 + const bool& etrue) const
  367 +{
  368 +
  369 + // Get effective area value in cm2
  370 + double aeff_max = m_logE.interpolate(logE, m_aeff);
  371 +
  372 + // Make sure that effective area is not negative
  373 + if (aeff_max < 0.0) {
  374 + aeff_max = 0.0;
  375 + }
  376 +
  377 + // Return result
  378 + return aeff_max;
  379 +
  380 +}
  381 +
354 382 /*==========================================================================
355 383 = =
356 384 = Private methods =
... ...
inst/cta/src/GCTAModelAeffBackground.cpp 0 → 100644
  View file @caa5b1e
  1 +/***************************************************************************
  2 + * GCTAModelAeffBackground.cpp - CTA Aeff background model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 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 GCTAModelAeffBackground.cpp
  23 + * @brief CTA Aeff background model class implementation
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GTools.hpp"
  32 +#include "GMath.hpp"
  33 +#include "GIntegral.hpp"
  34 +#include "GModelRegistry.hpp"
  35 +#include "GModelSpectralRegistry.hpp"
  36 +#include "GModelTemporalRegistry.hpp"
  37 +#include "GModelTemporalConst.hpp"
  38 +#include "GModelSpectralNodes.hpp"
  39 +#include "GCTAModelAeffBackground.hpp"
  40 +#include "GCTAObservation.hpp"
  41 +#include "GCTAResponseIrf.hpp"
  42 +#include "GCTABackground.hpp"
  43 +
  44 +/* __ Globals ____________________________________________________________ */
  45 +const GCTAModelAeffBackground g_cta_aeff_background_seed;
  46 +const GModelRegistry g_cta_aeff_background_registry(&g_cta_aeff_background_seed);
  47 +
  48 +/* __ Method name definitions ____________________________________________ */
  49 +#define G_EVAL "GCTAModelAeffBackground::eval(GEvent&, GObservation&)"
  50 +#define G_EVAL_GRADIENTS "GCTAModelAeffBackground::eval_gradients(GEvent&,"\
  51 + " GObservation&)"
  52 +#define G_NPRED "GCTAModelAeffBackground::npred(GEnergy&, GTime&,"\
  53 + " GObservation&)"
  54 +#define G_MC "GCTAModelAeffBackground::mc(GObservation&, GRan&)"
  55 +#define G_XML_SPECTRAL "GCTAModelAeffBackground::xml_spectral(GXmlElement&)"
  56 +#define G_XML_TEMPORAL "GCTAModelAeffBackground::xml_temporal(GXmlElement&)"
  57 +
  58 +/* __ Macros _____________________________________________________________ */
  59 +
  60 +/* __ Coding definitions _________________________________________________ */
  61 +#define G_USE_NPRED_CACHE
  62 +
  63 +/* __ Debug definitions __________________________________________________ */
  64 +//#define G_DUMP_MC
  65 +//#define G_DEBUG_NPRED
  66 +
  67 +/* __ Constants __________________________________________________________ */
  68 +
  69 +
  70 +/*==========================================================================
  71 + = =
  72 + = Constructors/destructors =
  73 + = =
  74 + ==========================================================================*/
  75 +
  76 +/***********************************************************************//**
  77 + * @brief Void constructor
  78 + ***************************************************************************/
  79 +GCTAModelAeffBackground::GCTAModelAeffBackground(void) : GModelData()
  80 +{
  81 + // Initialise class members
  82 + init_members();
  83 +
  84 + // Return
  85 + return;
  86 +}
  87 +
  88 +
  89 +/***********************************************************************//**
  90 + * @brief Constructor
  91 + *
  92 + * @param[in] xml XML element.
  93 + *
  94 + * Constructs a CTA effective area background model from the information
  95 + * provided by an XML element. The XML element is expected to have the
  96 + * following structure
  97 + *
  98 + * <source name="..." type="..." instrument="...">
  99 + * <spectrum type="...">
  100 + * ...
  101 + * </spectrum>
  102 + * </source>
  103 + *
  104 + * Optionally, a temporal model may be provided using the following
  105 + * syntax
  106 + *
  107 + * <source name="..." type="..." instrument="...">
  108 + * <spectrum type="...">
  109 + * ...
  110 + * </spectrum>
  111 + * <temporalModel type="...">
  112 + * ...
  113 + * </temporalModel>
  114 + * </source>
  115 + *
  116 + * If no temporal component is found a constant model is assumed.
  117 + ***************************************************************************/
  118 +GCTAModelAeffBackground::GCTAModelAeffBackground(const GXmlElement& xml) :
  119 + GModelData(xml)
  120 +{
  121 + // Initialise members
  122 + init_members();
  123 +
  124 + // Read XML
  125 + read(xml);
  126 +
  127 + // Set parameter pointers
  128 + set_pointers();
  129 +
  130 + // Return
  131 + return;
  132 +}
  133 +
  134 +
  135 +/***********************************************************************//**
  136 + * @brief Construct from spectral component
  137 + *
  138 + * @param[in] spectral Spectral model component.
  139 + *
  140 + * Constructs a CTA effective area background model from a spectral
  141 + * model component. The temporal component is assumed to be constant.
  142 + * Please refer to the class GModelSpectral to learn more about the
  143 + * definition of the spectral components.
  144 + ***************************************************************************/
  145 +GCTAModelAeffBackground::GCTAModelAeffBackground(const GModelSpectral& spectral) :
  146 + GModelData()
  147 +{
  148 + // Initialise members
  149 + init_members();
  150 +
  151 + // Allocate temporal constant model
  152 + GModelTemporalConst temporal;
  153 +
  154 + // Clone model components
  155 + m_spectral = spectral.clone();
  156 + m_temporal = temporal.clone();
  157 +
  158 + // Set parameter pointers
  159 + set_pointers();
  160 +
  161 + // Return
  162 + return;
  163 +}
  164 +
  165 +
  166 +/***********************************************************************//**
  167 + * @brief Copy constructor
  168 + *
  169 + * @param[in] aeff CTA instrument background model.
  170 + ***************************************************************************/
  171 +GCTAModelAeffBackground::GCTAModelAeffBackground(const GCTAModelAeffBackground& aeff) :
  172 + GModelData(aeff)
  173 +{
  174 + // Initialise class members
  175 + init_members();
  176 +
  177 + // Copy members
  178 + copy_members(aeff);
  179 +
  180 + // Return
  181 + return;
  182 +}
  183 +
  184 +
  185 +/***********************************************************************//**
  186 + * @brief Destructor
  187 + ***************************************************************************/
  188 +GCTAModelAeffBackground::~GCTAModelAeffBackground(void)
  189 +{
  190 + // Free members
  191 + free_members();
  192 +
  193 + // Return
  194 + return;
  195 +}
  196 +
  197 +
  198 +/*==========================================================================
  199 + = =
  200 + = Operators =
  201 + = =
  202 + ==========================================================================*/
  203 +
  204 +/***********************************************************************//**
  205 + * @brief Assignment operator
  206 + *
  207 + * @param[in] aeff CTA effective area background model.
  208 + * @return CTA effective area background model.
  209 + ***************************************************************************/
  210 +GCTAModelAeffBackground& GCTAModelAeffBackground::operator=(const GCTAModelAeffBackground& aeff)
  211 +{
  212 + // Execute only if object is not identical
  213 + if (this != &aeff) {
  214 +
  215 + // Copy base class members
  216 + this->GModelData::operator=(aeff);
  217 +
  218 + // Free members
  219 + free_members();
  220 +
  221 + // Initialise private members
  222 + init_members();
  223 +
  224 + // Copy members
  225 + copy_members(aeff);
  226 +
  227 + } // endif: object was not identical
  228 +
  229 + // Return this object
  230 + return *this;
  231 +}
  232 +
  233 +
  234 +/*==========================================================================
  235 + = =
  236 + = Public methods =
  237 + = =
  238 + ==========================================================================*/
  239 +
  240 +/***********************************************************************//**
  241 + * @brief Clear CTA effective area background model
  242 + *
  243 + * This method properly resets the CTA effective area background model to an
  244 + * initial state.
  245 + ***************************************************************************/
  246 +void GCTAModelAeffBackground::clear(void)
  247 +{
  248 + // Free class members (base and derived classes, derived class first)
  249 + free_members();
  250 + this->GModelData::free_members();
  251 +
  252 + // Initialise members
  253 + this->GModelData::init_members();
  254 + init_members();
  255 +
  256 + // Return
  257 + return;
  258 +}
  259 +
  260 +
  261 +/***********************************************************************//**
  262 + * @brief Clone CTA effective area background model
  263 + *
  264 + * @return Pointer to deep copy of CTA effective area background model.
  265 + ***************************************************************************/
  266 +GCTAModelAeffBackground* GCTAModelAeffBackground::clone(void) const
  267 +{
  268 + return new GCTAModelAeffBackground(*this);
  269 +}
  270 +
  271 +
  272 +/***********************************************************************//**
  273 + * @brief Evaluate function
  274 + *
  275 + * @param[in] event Observed event.
  276 + * @param[in] obs Observation.
  277 + * @return Function value.
  278 + *
  279 + * @exception GException::invalid_argument
  280 + * Specified observation is not of the expected type.
  281 + *
  282 + * @todo Make sure that DETX and DETY are always set in GCTAInstDir.
  283 + ***************************************************************************/
  284 +double GCTAModelAeffBackground::eval(const GEvent& event,
  285 + const GObservation& obs) const
  286 +{
  287 + // Get pointer on CTA observation
  288 + const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
  289 + if (cta == NULL) {
  290 + std::string msg = "Specified observation is not a CTA observation.\n" +
  291 + obs.print();
  292 + throw GException::invalid_argument(G_EVAL, msg);
  293 + }
  294 +
  295 + // Get pointer on CTA IRF response
  296 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
  297 + if (rsp == NULL) {
  298 + std::string msg = "Specified observation does not contain an IRF response.\n" +
  299 + obs.print();
  300 + throw GException::invalid_argument(G_EVAL, msg);
  301 + }
  302 +
  303 + // Retrieve pointer to CTA Effective Area
  304 + const GCTAAeff* aeff = rsp->aeff();
  305 + if (aeff == NULL) {
  306 + std::string msg = "Specified observation contains no effective area"
  307 + " information.\n" + obs.print();
  308 + throw GException::invalid_argument(G_EVAL, msg);
  309 + }
  310 +
  311 + // Extract CTA instrument direction from event
  312 + const GCTAInstDir* dir = dynamic_cast<const GCTAInstDir*>(&(event.dir()));
  313 + if (dir == NULL) {
  314 + std::string msg = "No CTA instrument direction found in event.";
  315 + throw GException::invalid_argument(G_EVAL, msg);
  316 + }
  317 +
  318 + // Set theta and phi from instrument coordinates
  319 + double theta = std::sqrt(dir->detx() * dir->detx() + dir->dety() * dir->dety());
  320 + double phi = gammalib::atan2d(dir->dety(), dir->detx()) * gammalib::deg2rad;
  321 +
  322 + // Evaluate function
  323 + double logE = event.energy().log10TeV();
  324 + double spat = (*aeff)(logE, theta, phi, cta->pointing().zenith(), cta->pointing().azimuth(), false);
  325 + double spec = (spectral() != NULL)
  326 + ? spectral()->eval(event.energy(), event.time()) : 1.0;
  327 + double temp = (temporal() != NULL)
  328 + ? temporal()->eval(event.time()) : 1.0;
  329 +
  330 + // Compute value
  331 + double value = spat * spec * temp;
  332 +
  333 + // Apply deadtime correction
  334 + value *= obs.deadc(event.time());
  335 +
  336 + // Return value
  337 + return value;
  338 +}
  339 +
  340 +
  341 +/***********************************************************************//**
  342 + * @brief Evaluate function and gradients
  343 + *
  344 + * @param[in] event Observed event.
  345 + * @param[in] obs Observation.
  346 + * @return Function value.
  347 + *
  348 + * @exception GException::invalid_argument
  349 + * Specified observation is not of the expected type.
  350 + *
  351 + * @todo Make sure that DETX and DETY are always set in GCTAInstDir.
  352 + ***************************************************************************/
  353 +double GCTAModelAeffBackground::eval_gradients(const GEvent& event,
  354 + const GObservation& obs) const
  355 +{
  356 + // Get pointer on CTA observation
  357 + const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
  358 + if (cta == NULL) {
  359 + std::string msg = "Specified observation is not a CTA observation.\n" +
  360 + obs.print();
  361 + throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
  362 + }
  363 +
  364 + // Get pointer on CTA IRF response
  365 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
  366 + if (rsp == NULL) {
  367 + std::string msg = "Specified observation does not contain an"
  368 + " IRF response.\n" + obs.print();
  369 + throw GException::invalid_argument(G_EVAL_GRADIENTS, msg);
  370 + }
  371 +
  372 + // Retrieve pointer to CTA Effective Area
  373 + const GCTAAeff* aeff = rsp->aeff();
  374 + if (aeff == NULL) {
  375 + std::string msg = "Specified observation contains no effective area"
  376 + " information.\n" + obs.print();
  377 + throw GException::invalid_argument(G_EVAL, msg);
  378 + }
  379 +
  380 + // Extract CTA instrument direction from event
  381 + const GCTAInstDir* dir = dynamic_cast<const GCTAInstDir*>(&(event.dir()));
  382 + if (dir == NULL) {
  383 + std::string msg = "No CTA instrument direction found in event.";
  384 + throw GException::invalid_argument(G_EVAL, msg);
  385 + }
  386 +
  387 + // Set theta and phi from instrument coordinates
  388 + double theta = std::sqrt(dir->detx() * dir->detx() + dir->dety() * dir->dety());
  389 + double phi = gammalib::atan2d(dir->dety() * gammalib::rad2deg, dir->detx() * gammalib::rad2deg) * gammalib::deg2rad;
  390 +
  391 + // Evaluate function
  392 + double logE = event.energy().log10TeV();
  393 + double spat = (*aeff)(logE, theta, phi, cta->pointing().zenith(), cta->pointing().azimuth(), false);
  394 + double spec = (spectral() != NULL)
  395 + ? spectral()->eval(event.energy(), event.time()) : 1.0;
  396 + double temp = (temporal() != NULL)
  397 + ? temporal()->eval(event.time()) : 1.0;
  398 +
  399 + // Compute value
  400 + double value = spat * spec * temp;
  401 +
  402 + // Apply deadtime correction
  403 + double deadc = obs.deadc(event.time());
  404 + value *= deadc;
  405 +
  406 + // Multiply factors to spectral gradients
  407 + if (spectral() != NULL) {
  408 + double fact = spat * temp * deadc;
  409 + if (fact != 1.0) {
  410 + for (int i = 0; i < spectral()->size(); ++i)
  411 + (*spectral())[i].factor_gradient( (*spectral())[i].factor_gradient() * fact );
  412 + }
  413 + }
  414 +
  415 + // Multiply factors to temporal gradients
  416 + if (temporal() != NULL) {
  417 + double fact = spat * spec * deadc;
  418 + if (fact != 1.0) {
  419 + for (int i = 0; i < temporal()->size(); ++i)
  420 + (*temporal())[i].factor_gradient( (*temporal())[i].factor_gradient() * fact );
  421 + }
  422 + }
  423 +
  424 + // Return value
  425 + return value;
  426 +}
  427 +
  428 +
  429 +/***********************************************************************//**
  430 + * @brief Return spatially integrated background model
  431 + *
  432 + * @param[in] obsEng Measured event energy.
  433 + * @param[in] obsTime Measured event time.
  434 + * @param[in] obs Observation.
  435 + * @return Spatially integrated model.
  436 + *
  437 + * @exception GException::invalid_argument
  438 + * The specified observation is not a CTA observation.
  439 + *
  440 + * Spatially integrates the instrumental background model for a given
  441 + * measured event energy and event time. This method also applies a deadtime
  442 + * correction factor, so that the normalization of the model is a real rate
  443 + * (counts/MeV/s).
  444 + ***************************************************************************/
  445 +double GCTAModelAeffBackground::npred(const GEnergy& obsEng,
  446 + const GTime& obsTime,
  447 + const GObservation& obs) const
  448 +{
  449 + // Set number of iterations for Romberg integration.
  450 + static const int iter_theta = 6;
  451 + static const int iter_phi = 6;
  452 +
  453 + // Initialise result
  454 + double npred = 0.0;
  455 + bool has_npred = false;
  456 +
  457 + // Build unique identifier
  458 + std::string id = obs.instrument() + "::" + obs.id();
  459 +
  460 + // Check if Npred value is already in cache
  461 + #if defined(G_USE_NPRED_CACHE)
  462 + if (!m_npred_names.empty()) {
  463 +
  464 + // Search for unique identifier, and if found, recover Npred value
  465 + // and break
  466 + for (int i = 0; i < m_npred_names.size(); ++i) {
  467 + if (m_npred_names[i] == id && m_npred_energies[i] == obsEng) {
  468 + npred = m_npred_values[i];
  469 + has_npred = true;
  470 + #if defined(G_DEBUG_NPRED)
  471 + std::cout << "GCTAModelAeffBackground::npred:";
  472 + std::cout << " cache=" << i;
  473 + std::cout << " npred=" << npred << std::endl;
  474 + #endif
  475 + break;
  476 + }
  477 + }
  478 +
  479 + } // endif: there were values in the Npred cache
  480 + #endif
  481 +
  482 + // Continue only if no Npred cache value has been found
  483 + if (!has_npred) {
  484 +
  485 + // Evaluate only if model is valid
  486 + if (valid_model()) {
  487 +
  488 + // Get log10 of energy in TeV
  489 + double logE = obsEng.log10TeV();
  490 +
  491 + // Spatially integrate radial component
  492 + npred = this->aeff_integral(obs, logE);
  493 +
  494 + // Store result in Npred cache
  495 + #if defined(G_USE_NPRED_CACHE)
  496 + m_npred_names.push_back(id);
  497 + m_npred_energies.push_back(obsEng);
  498 + m_npred_times.push_back(obsTime);
  499 + m_npred_values.push_back(npred);
  500 + #endif
  501 +
  502 + // Debug: Check for NaN
  503 + #if defined(G_NAN_CHECK)
  504 + if (gammalib::is_notanumber(npred) || gammalib::is_infinite(npred)) {
  505 + std::string origin = "GCTAModelAeffBackground::npred";
  506 + std::string message = " NaN/Inf encountered (npred=" +
  507 + gammalib::str(npred) + ")";
  508 + gammalib::warning(origin, message);
  509 + }
  510 + #endif
  511 +
  512 + } // endif: model was valid
  513 +
  514 + } // endif: Npred computation required
  515 +
  516 + // Multiply in spectral and temporal components
  517 + npred *= spectral()->eval(obsEng, obsTime);
  518 + npred *= temporal()->eval(obsTime);
  519 +
  520 + // Apply deadtime correction
  521 + npred *= obs.deadc(obsTime);
  522 +
  523 + // Return Npred
  524 + return npred;
  525 +}
  526 +
  527 +
  528 +/***********************************************************************//**
  529 + * @brief Return simulated list of events
  530 + *
  531 + * @param[in] obs Observation.
  532 + * @param[in] ran Random number generator.
  533 + * @return Pointer to list of simulated events (needs to be de-allocated by
  534 + * client)
  535 + *
  536 + * @exception GException::invalid_argument
  537 + * Specified observation is not a CTA observation.
  538 + *
  539 + * Draws a sample of events from the background model using a Monte
  540 + * Carlo simulation. The region of interest, the energy boundaries and the
  541 + * good time interval for the sampling will be extracted from the observation
  542 + * argument that is passed to the method. The method also requires a random
  543 + * number generator of type GRan which is passed by reference, hence the
  544 + * state of the random number generator will be changed by the method.
  545 + *
  546 + * The method also applies a deadtime correction using a Monte Carlo process,
  547 + * taking into account temporal deadtime variations. For this purpose, the
  548 + * method makes use of the time dependent GObservation::deadc method.
  549 + *
  550 + * For each event in the returned event list, the sky direction, the nominal
  551 + * coordinates (DETX and DETY), the energy and the time will be set.
  552 + ***************************************************************************/
  553 +GCTAEventList* GCTAModelAeffBackground::mc(const GObservation& obs, GRan& ran) const
  554 +{
  555 + // Initialise new event list
  556 + GCTAEventList* list = new GCTAEventList;
  557 +
  558 + // Continue only if model is valid)
  559 + if (valid_model()) {
  560 +
  561 + // Retrieve CTA observation
  562 + const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
  563 + if (cta == NULL) {
  564 + std::string msg = "Specified observation is not a CTA observation.\n" +
  565 + obs.print();
  566 + throw GException::invalid_argument(G_MC, msg);
  567 + }
  568 +
  569 + // Get pointer on CTA IRF response
  570 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
  571 + if (rsp == NULL) {
  572 + std::string msg = "Specified observation does not contain"
  573 + " an IRF response.\n" + obs.print();
  574 + throw GException::invalid_argument(G_MC, msg);
  575 + }
  576 +
  577 + // Retrieve CTA response and pointing
  578 + const GCTAPointing& pnt = cta->pointing();
  579 +
  580 + // Get pointer to CTA effective area
  581 + const GCTAAeff* aeff = rsp->aeff();
  582 + if (aeff == NULL) {
  583 + std::string msg = "Specified observation contains no effective area"
  584 + " information.\n" + obs.print();
  585 + throw GException::invalid_argument(G_MC, msg);
  586 + }
  587 +
  588 + // Retrieve event list to access the ROI, energy boundaries and GTIs
  589 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  590 + if (events == NULL) {
  591 + std::string msg = "No CTA event list found in observation.\n" +
  592 + obs.print();
  593 + throw GException::invalid_argument(G_MC, msg);
  594 + }
  595 +
  596 + // Get simulation region
  597 + const GCTARoi& roi = events->roi();
  598 + const GEbounds& ebounds = events->ebounds();
  599 + const GGti& gti = events->gti();
  600 +
  601 + // Get maximum offset value for simulations
  602 + double max_theta = pnt.dir().dist(roi.centre().dir()) + roi.radius() * gammalib::deg2rad;
  603 +
  604 + // Set simulation region for result event list
  605 + list->roi(roi);
  606 + list->ebounds(ebounds);
  607 + list->gti(gti);
  608 +
  609 + GEbounds spectral_ebounds = GEbounds(20, ebounds.emin(), ebounds.emax(), true);
  610 + GModelSpectralNodes spectral;
  611 + for (int i = 0; i < spectral_ebounds.size(); ++i) {
  612 + GEnergy energy = spectral_ebounds.elogmean(i);
  613 + double intensity = this->aeff_integral(obs, energy.log10TeV());
  614 + double norm = m_spectral->eval(energy, events->tstart());
  615 + spectral.append(energy, norm * intensity);
  616 + }
  617 +
  618 + // Loop over all energy boundaries
  619 + for (int ieng = 0; ieng < ebounds.size(); ++ieng) {
  620 +
  621 + // Compute the background rate in model within the energy
  622 + // boundaries from spectral component (units: cts/s).
  623 + // Note that the time here is ontime. Deadtime correction will
  624 + // be done later.
  625 + double rate = spectral.flux(ebounds.emin(ieng), ebounds.emax(ieng));
  626 +
  627 + // Debug option: dump rate
  628 + #if defined(G_DUMP_MC)
  629 + std::cout << "GCTAModelAeffBackground::mc(\"" << name() << "\": ";
  630 + std::cout << "rate=" << rate << " cts/s)" << std::endl;
  631 + #endif
  632 +
  633 + // Loop over all good time intervals
  634 + for (int itime = 0; itime < gti.size(); ++itime) {
  635 +
  636 + // Get Monte Carlo event arrival times from temporal model
  637 + GTimes times = m_temporal->mc(rate,
  638 + gti.tstart(itime),
  639 + gti.tstop(itime),
  640 + ran);
  641 +
  642 + // Get number of events
  643 + int n_events = times.size();
  644 +
  645 + // Reserve space for events
  646 + if (n_events > 0) {
  647 + list->reserve(n_events);
  648 + }
  649 +
  650 + // Debug option: provide number of times and initialize
  651 + // statisics
  652 + #if defined(G_DUMP_MC)
  653 + std::cout << " Interval " << itime;
  654 + std::cout << " times=" << n_events << std::endl;
  655 + int n_killed_by_deadtime = 0;
  656 + int n_killed_by_roi = 0;
  657 + #endif
  658 +
  659 + // Loop over events
  660 + for (int i = 0; i < n_events; ++i) {
  661 +
  662 + // Apply deadtime correction
  663 + double deadc = obs.deadc(times[i]);
  664 + if (deadc < 1.0) {
  665 + if (ran.uniform() > deadc) {
  666 + #if defined(G_DUMP_MC)
  667 + n_killed_by_deadtime++;
  668 + #endif
  669 + continue;
  670 + }
  671 + }
  672 +
  673 + // Get Monte Carlo event energy from spectral model
  674 + GEnergy energy = spectral.mc(ebounds.emin(ieng),
  675 + ebounds.emax(ieng),
  676 + times[i],
  677 + ran);
  678 +
  679 + // Get maximum effective area for rejection method
  680 + double max_aeff = aeff->max(energy.log10TeV(), pnt.zenith(), pnt.azimuth(), false);
  681 +
  682 + // Initialise randomised coordinates
  683 + double offset = 0.0;
  684 + double phi = 0.0;
  685 +
  686 + // Initialise acceptance fraction for rejection method
  687 + double acceptance_fraction = 0.0;
  688 +
  689 + // Start rejection method loop
  690 + do {
  691 +
  692 + // Throw random offset and azimuth angle in considered range
  693 + offset = std::sqrt(ran.uniform()) * max_theta;
  694 + phi = ran.uniform() * gammalib::twopi;
  695 +
  696 + // Compute function value at this offset angle
  697 + double value = (*aeff)(energy.log10TeV(), offset, phi, pnt.zenith(), pnt.azimuth(), false);
  698 +
  699 + // Compute acceptance fraction
  700 + acceptance_fraction = value / max_aeff;
  701 +
  702 + } while (ran.uniform() > acceptance_fraction);
  703 +
  704 +
  705 + // Convert CTA pointing direction in instrument system
  706 + GCTAInstDir mc_dir(pnt.dir());
  707 +
  708 + // Rotate pointing direction by offset and azimuth angle
  709 + mc_dir.dir().rotate_deg(phi * gammalib::rad2deg, offset * gammalib::rad2deg);
  710 +
  711 + // Compute DETX and DETY coordinates
  712 + double detx(0.0);
  713 + double dety(0.0);
  714 + if (offset > 0.0 ) {
  715 + detx = offset * std::cos(phi);
  716 + dety = offset * std::sin(phi);
  717 + }
  718 +
  719 + // Set DETX and DETY coordinates
  720 + mc_dir.detx(detx);
  721 + mc_dir.dety(dety);
  722 +
  723 + // Allocate event
  724 + GCTAEventAtom event;
  725 +
  726 + // Set event attributes
  727 + event.dir(mc_dir);
  728 + event.energy(energy);
  729 + event.time(times[i]);
  730 +
  731 + // Append event to list if it falls in ROI
  732 + if (events->roi().contains(event)) {
  733 + list->append(event);
  734 + }
  735 + #if defined(G_DUMP_MC)
  736 + else {
  737 + n_killed_by_roi++;
  738 + }
  739 + #endif
  740 +
  741 + } // endfor: looped over all events
  742 +
  743 + // Debug option: provide statisics
  744 + #if defined(G_DUMP_MC)
  745 + std::cout << " Killed by deadtime=";
  746 + std::cout << n_killed_by_deadtime << std::endl;
  747 + std::cout << " Killed by ROI=";
  748 + std::cout << n_killed_by_roi << std::endl;
  749 + #endif
  750 +
  751 + } // endfor: looped over all GTIs
  752 +
  753 + } // endfor: looped over all energy boundaries
  754 +
  755 + } // endif: model was valid
  756 +
  757 + // Return
  758 + return list;
  759 +}
  760 +
  761 +
  762 +/***********************************************************************//**
  763 + * @brief Read CTA instrument background model from XML element
  764 + *
  765 + * @param[in] xml XML element.
  766 + *
  767 + * Set up CTA instrument background model from the information provided by
  768 + * an XML element. The XML element is expected to have the following
  769 + * structure
  770 + *
  771 + * <source name="..." type="..." instrument="...">
  772 + * <spectrum type="...">
  773 + * ...
  774 + * </spectrum>
  775 + * </source>
  776 + *
  777 + * Optionally, a temporal model may be provided using the following
  778 + * syntax
  779 + *
  780 + * <source name="..." type="..." instrument="...">
  781 + * <spectrum type="...">
  782 + * ...
  783 + * </spectrum>
  784 + * <temporalModel type="...">
  785 + * ...
  786 + * </temporalModel>
  787 + * </source>
  788 + *
  789 + * If no temporal component is found a constant model is assumed.
  790 + ***************************************************************************/
  791 +void GCTAModelAeffBackground::read(const GXmlElement& xml)
  792 +{
  793 + // Clear model
  794 + clear();
  795 +
  796 + // Initialise XML elements
  797 + const GXmlElement* spectral = NULL;
  798 + const GXmlElement* temporal = NULL;
  799 +
  800 + // Get pointer on spectrum
  801 + spectral = xml.element("spectrum", 0);
  802 +
  803 + // Extract spectral model
  804 + m_spectral = xml_spectral(*spectral);
  805 +
  806 + // Optionally get temporal model
  807 + if (xml.elements("temporalModel")) {
  808 + temporal = xml.element("temporalModel", 0);
  809 + m_temporal = xml_temporal(*temporal);
  810 + }
  811 + else {
  812 + GModelTemporalConst constant;
  813 + m_temporal = constant.clone();
  814 + }
  815 +
  816 + // Set model name
  817 + name(xml.attribute("name"));
  818 +
  819 + // Set instruments
  820 + instruments(xml.attribute("instrument"));
  821 +
  822 + // Set observation identifiers
  823 + ids(xml.attribute("id"));
  824 +
  825 + // Check flag if TS value should be computed
  826 + bool tscalc = (xml.attribute("tscalc") == "1") ? true : false;
  827 +
  828 + // Set flag if TS value should be computed
  829 + this->tscalc(tscalc);
  830 +
  831 + // Set parameter pointers
  832 + set_pointers();
  833 +
  834 + // Return
  835 + return;
  836 +}
  837 +
  838 +
  839 +/***********************************************************************//**
  840 + * @brief Write CTA instrument background model into XML element
  841 + *
  842 + * @param[in] xml XML element.
  843 + *
  844 + * Write CTA instrument background model information into an XML element.
  845 + * The XML element will have the following structure
  846 + *
  847 + * <source name="..." type="..." instrument="...">
  848 + * <spectrum type="...">
  849 + * ...
  850 + * </spectrum>
  851 + * </source>
  852 + *
  853 + * If the model contains a non-constant temporal model, the temporal
  854 + * component will also be written following the syntax
  855 + *
  856 + * <source name="..." type="..." instrument="...">
  857 + * <spectrum type="...">
  858 + * ...
  859 + * </spectrum>
  860 + * <temporalModel type="...">
  861 + * ...
  862 + * </temporalModel>
  863 + * </source>
  864 + *
  865 + * If no temporal component is found a constant model is assumed.
  866 + ***************************************************************************/
  867 +void GCTAModelAeffBackground::write(GXmlElement& xml) const
  868 +{
  869 + // Initialise pointer on source
  870 + GXmlElement* src = NULL;
  871 +
  872 + // Search corresponding source
  873 + int n = xml.elements("source");
  874 + for (int k = 0; k < n; ++k) {
  875 + GXmlElement* element = xml.element("source", k);
  876 + if (element->attribute("name") == name()) {
  877 + src = element;
  878 + break;
  879 + }
  880 + }
  881 +
  882 + // If we have a temporal model that is either not a constant, or a
  883 + // constant with a normalization value that differs from 1.0 then
  884 + // write the temporal component into the XML element. This logic
  885 + // assures compatibility with the Fermi/LAT format as this format
  886 + // does not handle temporal components.
  887 + bool write_temporal = ((m_temporal != NULL) &&
  888 + (m_temporal->type() != "Constant" ||
  889 + (*m_temporal)[0].value() != 1.0));
  890 +
  891 + // If no source with corresponding name was found then append one
  892 + if (src == NULL) {
  893 + src = xml.append("source");
  894 + if (spectral() != NULL) src->append(GXmlElement("spectrum"));
  895 + if (write_temporal) src->append(GXmlElement("temporalModel"));
  896 + }
  897 +
  898 + // Set model type, name and optionally instruments
  899 + src->attribute("name", name());
  900 + src->attribute("type", type());
  901 + if (instruments().length() > 0) {
  902 + src->attribute("instrument", instruments());
  903 + }
  904 + std::string identifiers = ids();
  905 + if (identifiers.length() > 0) {
  906 + src->attribute("id", identifiers);
  907 + }
  908 +
  909 + // Write spectral model
  910 + if (spectral() != NULL) {
  911 + GXmlElement* spec = src->element("spectrum", 0);
  912 + spectral()->write(*spec);
  913 + }
  914 +
  915 + // Optionally write temporal model
  916 + if (write_temporal) {
  917 + if (dynamic_cast<GModelTemporalConst*>(temporal()) == NULL) {
  918 + GXmlElement* temp = src->element("temporalModel", 0);
  919 + temporal()->write(*temp);
  920 + }
  921 + }
  922 +
  923 + // Return
  924 + return;
  925 +}
  926 +
  927 +
  928 +/***********************************************************************//**
  929 + * @brief Print CTA instrument background model information
  930 + *
  931 + * @param[in] chatter Chattiness (defaults to NORMAL).
  932 + * @return String containing CTA instrument background model information.
  933 + ***************************************************************************/
  934 +std::string GCTAModelAeffBackground::print(const GChatter& chatter) const
  935 +{
  936 + // Initialise result string
  937 + std::string result;
  938 +
  939 + // Continue only if chatter is not silent
  940 + if (chatter != SILENT) {
  941 +
  942 + // Append header
  943 + result.append("=== GCTAModelAeffBackground ===");
  944 +
  945 + // Determine number of parameters per type
  946 + int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
  947 + int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
  948 +
  949 + // Append attributes
  950 + result.append("\n"+print_attributes());
  951 +
  952 + // Append model type
  953 + result.append("\n"+gammalib::parformat("Model type"));
  954 + if (n_spectral > 0) {
  955 + result.append("\""+spectral()->type()+"\"");
  956 + if (n_temporal > 0) {
  957 + result.append(" * ");
  958 + }
  959 + }
  960 + if (n_temporal > 0) {
  961 + result.append("\""+temporal()->type()+"\"");
  962 + }
  963 +
  964 + // Append parameters
  965 + result.append("\n"+gammalib::parformat("Number of parameters") +
  966 + gammalib::str(size()));
  967 + result.append("\n"+gammalib::parformat("Number of spectral par's") +
  968 + gammalib::str(n_spectral));
  969 + for (int i = 0; i < n_spectral; ++i) {
  970 + result.append("\n"+(*spectral())[i].print());
  971 + }
  972 + result.append("\n"+gammalib::parformat("Number of temporal par's") +
  973 + gammalib::str(n_temporal));
  974 + for (int i = 0; i < n_temporal; ++i) {
  975 + result.append("\n"+(*temporal())[i].print());
  976 + }
  977 +
  978 + } // endif: chatter was not silent
  979 +
  980 + // Return result
  981 + return result;
  982 +}
  983 +
  984 +
  985 +/*==========================================================================
  986 + = =
  987 + = Private methods =
  988 + = =
  989 + ==========================================================================*/
  990 +
  991 +/***********************************************************************//**
  992 + * @brief Initialise class members
  993 + ***************************************************************************/
  994 +void GCTAModelAeffBackground::init_members(void)
  995 +{
  996 + // Initialise members
  997 + m_spectral = NULL;
  998 + m_temporal = NULL;
  999 +
  1000 + // Initialise Npred cache
  1001 + m_npred_names.clear();
  1002 + m_npred_energies.clear();
  1003 + m_npred_times.clear();
  1004 + m_npred_values.clear();
  1005 +
  1006 + // Return
  1007 + return;
  1008 +}
  1009 +
  1010 +
  1011 +/***********************************************************************//**
  1012 + * @brief Copy class members
  1013 + *
  1014 + * @param[in] bgd CTA background model.
  1015 + ***************************************************************************/
  1016 +void GCTAModelAeffBackground::copy_members(const GCTAModelAeffBackground& bgd)
  1017 +{
  1018 + // Copy cache
  1019 + m_npred_names = bgd.m_npred_names;
  1020 + m_npred_energies = bgd.m_npred_energies;
  1021 + m_npred_times = bgd.m_npred_times;
  1022 + m_npred_values = bgd.m_npred_values;
  1023 +
  1024 + // Clone spectral and temporal model components
  1025 + m_spectral = (bgd.m_spectral != NULL) ? bgd.m_spectral->clone() : NULL;
  1026 + m_temporal = (bgd.m_temporal != NULL) ? bgd.m_temporal->clone() : NULL;
  1027 +
  1028 + // Set parameter pointers
  1029 + set_pointers();
  1030 +
  1031 + // Return
  1032 + return;
  1033 +}
  1034 +
  1035 +
  1036 +/***********************************************************************//**
  1037 + * @brief Delete class members
  1038 + ***************************************************************************/
  1039 +void GCTAModelAeffBackground::free_members(void)
  1040 +{
  1041 + // Free memory
  1042 + if (m_spectral != NULL) delete m_spectral;
  1043 + if (m_temporal != NULL) delete m_temporal;
  1044 +
  1045 + // Signal free pointers
  1046 + m_spectral = NULL;
  1047 + m_temporal = NULL;
  1048 +
  1049 + // Return
  1050 + return;
  1051 +}
  1052 +
  1053 +
  1054 +/***********************************************************************//**
  1055 + * @brief Set pointers
  1056 + *
  1057 + * Set pointers to all model parameters. The pointers are stored in a vector
  1058 + * that is member of the GModelData base class.
  1059 + ***************************************************************************/
  1060 +void GCTAModelAeffBackground::set_pointers(void)
  1061 +{
  1062 + // Clear parameters
  1063 + m_pars.clear();
  1064 +
  1065 + // Determine the number of parameters
  1066 + int n_spectral = (spectral() != NULL) ? spectral()->size() : 0;
  1067 + int n_temporal = (temporal() != NULL) ? temporal()->size() : 0;
  1068 + int n_pars = n_spectral + n_temporal;
  1069 +
  1070 + // Continue only if there are parameters
  1071 + if (n_pars > 0) {
  1072 +
  1073 + // Gather spectral parameters
  1074 + for (int i = 0; i < n_spectral; ++i) {
  1075 + m_pars.push_back(&((*spectral())[i]));
  1076 + }
  1077 +
  1078 + // Gather temporal parameters
  1079 + for (int i = 0; i < n_temporal; ++i) {
  1080 + m_pars.push_back(&((*temporal())[i]));
  1081 + }
  1082 +
  1083 + }
  1084 +
  1085 + // Return
  1086 + return;
  1087 +}
  1088 +
  1089 +
  1090 +/***********************************************************************//**
  1091 + * @brief Verifies if model has all components
  1092 + *
  1093 + * Returns 'true' if models has a spectral and a temporal component.
  1094 + * Otherwise returns 'false'.
  1095 + ***************************************************************************/
  1096 +bool GCTAModelAeffBackground::valid_model(void) const
  1097 +{
  1098 + // Set result
  1099 + bool result = ((spectral() != NULL) && (temporal() != NULL));
  1100 +
  1101 + // Return result
  1102 + return result;
  1103 +}
  1104 +
  1105 +
  1106 +/***********************************************************************//**
  1107 + * @brief Construct spectral model from XML element
  1108 + *
  1109 + * @param[in] spectral XML element containing spectral model information.
  1110 + *
  1111 + * @exception GException::model_invalid_spectral
  1112 + * Invalid spectral model type encountered.
  1113 + *
  1114 + * Returns pointer to a spectral model that is defined in an XML element.
  1115 + ***************************************************************************/
  1116 +GModelSpectral* GCTAModelAeffBackground::xml_spectral(const GXmlElement& spectral) const
  1117 +{
  1118 + // Get spectral model type
  1119 + std::string type = spectral.attribute("type");
  1120 +
  1121 + // Get spectral model
  1122 + GModelSpectralRegistry registry;
  1123 + GModelSpectral* ptr = registry.alloc(type);
  1124 +
  1125 + // If model if valid then read model from XML file
  1126 + if (ptr != NULL) {
  1127 + ptr->read(spectral);
  1128 + }
  1129 +
  1130 + // ... otherwise throw an exception
  1131 + else {
  1132 + throw GException::model_invalid_spectral(G_XML_SPECTRAL, type);
  1133 + }
  1134 +
  1135 + // Return pointer
  1136 + return ptr;
  1137 +}
  1138 +
  1139 +
  1140 +/***********************************************************************//**
  1141 + * @brief Construct temporal model from XML element
  1142 + *
  1143 + * @param[in] temporal XML element containing temporal model information.
  1144 + *
  1145 + * @exception GException::model_invalid_temporal
  1146 + * Invalid temporal model type encountered.
  1147 + *
  1148 + * Returns pointer to a temporal model that is defined in an XML element.
  1149 + ***************************************************************************/
  1150 +GModelTemporal* GCTAModelAeffBackground::xml_temporal(const GXmlElement& temporal) const
  1151 +{
  1152 + // Get temporal model type
  1153 + std::string type = temporal.attribute("type");
  1154 +
  1155 + // Get temporal model
  1156 + GModelTemporalRegistry registry;
  1157 + GModelTemporal* ptr = registry.alloc(type);
  1158 +
  1159 + // If model if valid then read model from XML file
  1160 + if (ptr != NULL) {
  1161 + ptr->read(temporal);
  1162 + }
  1163 +
  1164 + // ... otherwise throw an exception
  1165 + else {
  1166 + throw GException::model_invalid_temporal(G_XML_TEMPORAL, type);
  1167 + }
  1168 +
  1169 + // Return pointer
  1170 + return ptr;
  1171 +}
  1172 +
  1173 +
  1174 +double GCTAModelAeffBackground::aeff_integral(const GObservation& obs, const double& logE) const
  1175 +{
  1176 + // Initialise result
  1177 + double value = 0.0;
  1178 +
  1179 + // Set number of iterations for Romberg integration.
  1180 + static const int iter_theta = 6;
  1181 + static const int iter_phi = 6;
  1182 +
  1183 + // Get pointer on CTA observation
  1184 + const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(&obs);
  1185 + if (cta == NULL) {
  1186 + std::string msg = "Specified observation is not a CTA"
  1187 + " observation.\n" + obs.print();
  1188 + throw GException::invalid_argument(G_NPRED, msg);
  1189 + }
  1190 +
  1191 + // Get pointer on CTA IRF response
  1192 + const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
  1193 + if (rsp == NULL) {
  1194 + std::string msg = "Specified observation does not contain"
  1195 + " an IRF response.\n" + obs.print();
  1196 + throw GException::invalid_argument(G_NPRED, msg);
  1197 + }
  1198 +
  1199 + // Retrieve pointer to CTA effective area
  1200 + const GCTAAeff* aeff = rsp->aeff();
  1201 + if (aeff == NULL) {
  1202 + std::string msg = "Specified observation contains no effective area"
  1203 + " information.\n" + obs.print();
  1204 + throw GException::invalid_argument(G_NPRED, msg);
  1205 + }
  1206 +
  1207 + // Get CTA event list
  1208 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(obs.events());
  1209 + if (events == NULL) {
  1210 + std::string msg = "No CTA event list found in observation.\n" +
  1211 + obs.print();
  1212 + throw GException::invalid_argument(G_NPRED, msg);
  1213 + }
  1214 +
  1215 + // Get reference to ROI centre
  1216 + const GSkyDir& roi_centre = events->roi().centre().dir();
  1217 +
  1218 + // Get ROI radius in radians
  1219 + double roi_radius = events->roi().radius() * gammalib::deg2rad;
  1220 +
  1221 + // Setup integration function
  1222 + GCTAModelAeffBackground::npred_roi_kern_theta integrand(aeff,
  1223 + logE,
  1224 + iter_phi);
  1225 +
  1226 + // Setup integration
  1227 + GIntegral integral(&integrand);
  1228 +
  1229 + // Set fixed number of iterations
  1230 + integral.fixed_iter(iter_theta);
  1231 +
  1232 + // Spatially integrate radial component
  1233 + value = integral.romberg(0.0, roi_radius);
  1234 +
  1235 + // Debug: Check for NaN
  1236 + #if defined(G_NAN_CHECK)
  1237 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  1238 + std::string origin = "GCTAModelAeffBackground::aeff_integral";
  1239 + std::string message = " NaN/Inf encountered (value=" +
  1240 + gammalib::str(value) + ", roi_radius=" +
  1241 + gammalib::str(roi_radius) + ")";
  1242 + gammalib::warning(origin, message);
  1243 + }
  1244 + #endif
  1245 +
  1246 + // Return
  1247 + return value;
  1248 +
  1249 +}
  1250 +
  1251 +
  1252 +/***********************************************************************//**
  1253 + * @brief Kernel for offset angle integration of background model
  1254 + *
  1255 + * @param[in] theta Offset angle from ROI centre (radians).
  1256 + *
  1257 + * Computes
  1258 + *
  1259 + * \f[
  1260 + * K(\rho | E, t) = \sin \theta \times
  1261 + * \int_{0}^{2\pi}
  1262 + * B(\theta,\phi | E, t) d\phi
  1263 + * \f]
  1264 + *
  1265 + * where \f$B(\theta,\phi | E, t)\f$ is the background model for a specific
  1266 + * observed energy \f$E\f$ and time \f$t\f$.
  1267 + ***************************************************************************/
  1268 +double GCTAModelAeffBackground::npred_roi_kern_theta::eval(const double& theta)
  1269 +{
  1270 + // Initialise value
  1271 + double value = 0.0;
  1272 +
  1273 + // Continue only if offset angle is positive
  1274 + if (theta > 0.0) {
  1275 +
  1276 + // Setup phi integration kernel
  1277 + GCTAModelAeffBackground::npred_roi_kern_phi integrand(m_aeff,
  1278 + m_logE,
  1279 + theta);
  1280 +
  1281 + // Setup integration
  1282 + GIntegral integral(&integrand);
  1283 +
  1284 + // Set fixed number of iterations
  1285 + integral.fixed_iter(m_iter);
  1286 +
  1287 + // Integrate over phi
  1288 + value = integral.romberg(0.0, gammalib::twopi) * std::sin(theta);
  1289 +
  1290 + // Debug: Check for NaN
  1291 + #if defined(G_NAN_CHECK)
  1292 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  1293 + std::string origin = "GCTAModelAeffBackground::npred_roi_kern_theta::eval"
  1294 + "(" + gammalib::str(theta) + ")";
  1295 + std::string message = " NaN/Inf encountered (value=" +
  1296 + gammalib::str(value) + ")";
  1297 + gammalib::warning(origin, message);
  1298 + }
  1299 + #endif
  1300 +
  1301 + } // endif: offset angle was positive
  1302 +
  1303 + // Return value
  1304 + return value;
  1305 +}
  1306 +
  1307 +
  1308 +/***********************************************************************//**
  1309 + * @brief Kernel for azimuth angle integration of background model
  1310 + *
  1311 + * @param[in] phi Azimuth angle around ROI centre (radians).
  1312 + *
  1313 + * Computes
  1314 + *
  1315 + * \f[
  1316 + * B(\theta, \phi | E, t)
  1317 + * \f]
  1318 + *
  1319 + * using
  1320 + *
  1321 + * \f[ {\rm detx} = \theta \cos \phi \f]
  1322 + * \f[ {\rm dety} = \theta \sin \phi \f]
  1323 + *
  1324 + * @todo Verify correct orientation of detx and dety with respect to phi
  1325 + ***************************************************************************/
  1326 +double GCTAModelAeffBackground::npred_roi_kern_phi::eval(const double& phi)
  1327 +{
  1328 + // Get background value
  1329 + double value = (*m_aeff)(m_logE, m_theta, phi);
  1330 +
  1331 + // Debug: Check for NaN
  1332 + #if defined(G_NAN_CHECK)
  1333 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  1334 + std::string origin = "GCTAModelAeffBackground::npred_roi_kern_phi::eval"
  1335 + "(" + gammalib::str(phi) + ")";
  1336 + std::string message = " NaN/Inf encountered (value=" +
  1337 + gammalib::str(value) + ", theta=" +
  1338 + gammalib::str(m_theta) + ", phi=" +
  1339 + gammalib::str(phi) + ")";
  1340 + gammalib::warning(origin, message);
  1341 + }
  1342 + #endif
  1343 +
  1344 + // Return Npred
  1345 + return value;
  1346 +}
... ...
inst/cta/test/data/cta_model_aeff_bgd.xml 0 → 100644
  View file @caa5b1e
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="My model" type="CTAAeffBackground" instrument="CTA">
  4 + <spectrum type="PowerLaw">
  5 + <parameter name="Prefactor" scale="1e-14" value="1.0" min="1e-3" max="1e3" free="1"/>
  6 + <parameter name="Index" scale="1.0" value="-2.7" min="-5.0" max="+5.0" free="1"/>
  7 + <parameter name="Scale" scale="1e6" value="1.0" min="0.01" max="1000.0" free="0"/>
  8 + </spectrum>
  9 + </source>
  10 +</source_library>
... ...
inst/cta/test/test_CTA.cpp
  View file @caa5b1e
... ... @@ -58,6 +58,7 @@ const std::string cta_model_xml = datadir+&quot;/crab.xml&quot;;
58 58 const std::string cta_rsp_xml = datadir+"/rsp_models.xml";
59 59 const std::string cta_cube_bgd_xml = datadir+"/cta_model_cube_bgd.xml";
60 60 const std::string cta_irf_bgd_xml = datadir+"/cta_model_irf_bgd.xml";
  61 +const std::string cta_aeff_bgd_xml = datadir+"/cta_model_aeff_bgd.xml";
61 62 const std::string cta_caldb_king = PACKAGE_SOURCE"/inst/cta/caldb/data/cta/e/bcf/IFAE20120510_50h_King";
62 63 const std::string cta_irf_king = "irf_file.fits";
63 64 const std::string cta_edisp_perf = PACKAGE_SOURCE"/inst/cta/test/caldb/cta_dummy_irf.dat";
... ... @@ -120,6 +121,7 @@ void TestGCTAModel::set(void)
120 121 // Append tests to test suite
121 122 append(static_cast<pfunction>(&TestGCTAModel::test_model_cube_bgd), "Test CTA cube background model");
122 123 append(static_cast<pfunction>(&TestGCTAModel::test_model_irf_bgd), "Test CTA IRF background model");
  124 + append(static_cast<pfunction>(&TestGCTAModel::test_model_aeff_bgd), "Test CTA Aeff background model");
123 125  
124 126 // Return
125 127 return;
... ... @@ -1198,6 +1200,74 @@ void TestGCTAModel::test_model_irf_bgd(void)
1198 1200  
1199 1201  
1200 1202 /***********************************************************************//**
  1203 + * @brief Test CTA IRF background model
  1204 + ***************************************************************************/
  1205 +void TestGCTAModel::test_model_aeff_bgd(void)
  1206 +{
  1207 + // Test void constuctor
  1208 + test_try("Test void constuctor");
  1209 + try {
  1210 + GCTAModelAeffBackground model;
  1211 + test_try_success();
  1212 + }
  1213 + catch (std::exception &e) {
  1214 + test_try_failure(e);
  1215 + }
  1216 +
  1217 + // Test XML constuctor
  1218 + test_try("Test XML constuctor");
  1219 + try {
  1220 + GXml xml(cta_aeff_bgd_xml);
  1221 + const GXmlElement& lib = *xml.element("source_library", 0);
  1222 + const GXmlElement& src = *lib.element("source", 0);
  1223 + GCTAModelIrfBackground model(src);
  1224 + test_value(model["Prefactor"].value(), 1.0);
  1225 + test_value(model["Index"].value(), 0.0);
  1226 + test_value(model["PivotEnergy"].value(), 1.0e6);
  1227 + test_assert(model.is_constant(), "Model is expected to be constant.");
  1228 + test_try_success();
  1229 + }
  1230 + catch (std::exception &e) {
  1231 + test_try_failure(e);
  1232 + }
  1233 +
  1234 + // Test spectral constuctor
  1235 + test_try("Test spectral constuctor");
  1236 + try {
  1237 + GModelSpectralPlaw plaw(1.0, 0.0, GEnergy(1.0, "TeV"));
  1238 + GCTAModelAeffBackground model(plaw);
  1239 + test_value(model["Prefactor"].value(), 1.0);
  1240 + test_value(model["Index"].value(), 0.0);
  1241 + test_value(model["PivotEnergy"].value(), 1.0e6);
  1242 + test_assert(model.is_constant(), "Model is expected to be constant.");
  1243 + test_try_success();
  1244 + }
  1245 + catch (std::exception &e) {
  1246 + test_try_failure(e);
  1247 + }
  1248 +
  1249 + // Test XML loading of instrumental background
  1250 + GModels models(cta_aeff_bgd_xml);
  1251 + GModel* model = models["My model"];
  1252 + test_value((*model)["Prefactor"].value(), 1.0);
  1253 + test_value((*model)["Index"].value(), 0.0);
  1254 + test_value((*model)["PivotEnergy"].value(), 1.0e6);
  1255 + test_assert(model->is_constant(), "Model is expected to be constant.");
  1256 +
  1257 + // Test XML saving and reloading of instrumental background
  1258 + models.save("test.xml");
  1259 + models.load("test.xml");
  1260 + model = models["My model"];
  1261 + test_value((*model)["Prefactor"].value(), 1.0);
  1262 + test_value((*model)["Index"].value(), 0.0);
  1263 + test_value((*model)["PivotEnergy"].value(), 1.0e6);
  1264 + test_assert(model->is_constant(), "Model is expected to be constant.");
  1265 +
  1266 + // Return
  1267 + return;
  1268 +}
  1269 +
  1270 +/***********************************************************************//**
1201 1271 * @brief Test unbinned observation handling
1202 1272 ***************************************************************************/
1203 1273 void TestGCTAObservation::test_unbinned_obs(void)
... ...
inst/cta/test/test_CTA.hpp
  View file @caa5b1e
... ... @@ -92,6 +92,7 @@ public:
92 92 virtual std::string classname(void) const { return "TestGCTAModel"; }
93 93 void test_model_cube_bgd(void);
94 94 void test_model_irf_bgd(void);
  95 + void test_model_aeff_bgd(void);
95 96 };
96 97  
97 98  
... ...