Sergio Hernández Cadena / gammalib

Commit 8e14f9ee5f050e07f21173d6c0c607fae05ebfd1

Authored by Jürgen Knödlseder
1 parent ad26baa0

Add GModelSpectralBins class (#3606)

Inline Side-by-side
Showing 16 changed files with 2057 additions and 13 deletions
ChangeLog
  View file @8e14f9e
1   -2021-04-02
  1 +2021-04-13
2 2  
3 3 * Version 2.0.0 released
4 4 ========================
5 5  
  6 + Add GModelSpectralBins class (#3606)
6 7 Add carbon footprint estimation to GApplication (#3592)
7 8 Add phase selection to GCOMSelection and GCOMDri (#3564)
8 9 Fix incompatibility in GMatrixSparse::mix_column_prepare() method (#3563)
... ...
1 1 New Features and Important Changes in GammaLib 2.0.0
2 2  
3   -2 April 2021
  3 +13 April 2021
4 4  
5 5  
6 6 1. Introduction
... ... @@ -19,6 +19,7 @@ The following classes have been added:
19 19 - GSkyRegionRect
20 20 - GModelAssociation
21 21 - GModelAssociations
  22 +- GModelSpectralBins
22 23 - GCTAModelSkyCube
23 24 - GCTABackground2D
24 25  
... ... @@ -159,8 +160,11 @@ Implemented GMatrixSparse::row() and GMatrixSparse::add_to_row() methods.
159 160  
160 161 8. Model module
161 162 ---------------
  163 +Added GModelSpectralBins class that implements statistically independent spectral
  164 +bins as the spectral model component (#3606).
  165 +
162 166 Added methods to GModelSpectralFunc to allow building of a file function
163   -in memort (#3535).
  167 +in memory (#3535).
164 168  
165 169 Added GModelAssociation class to implement a model association and the
166 170 GModelAssociations container class to hold a list of model associations.
... ...
README.md
  View file @8e14f9e
1 1 GammaLib information
2 2 ====================
3   -* Version: 2.0.0.dev (2 April 2021)
  3 +* Version: 2.0.0.dev (13 April 2021)
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  
... ... @@ -194,7 +194,7 @@ extensive unit test by typing:
194 194 If everything works successfully you should see
195 195  
196 196 ===================
197   - All 21 tests passed
  197 + All 22 tests passed
198 198 ===================
199 199  
200 200 or
... ...
doc/source/admin/release_history/2.0.rst
  View file @8e14f9e
1 1 .. _2.0:
2 2  
3   -GammaLib 2.0 release (13 March 2021)
  3 +GammaLib 2.0 release (13 April 2021)
4 4 ====================================
5 5  
6 6 Overview
... ... @@ -73,6 +73,8 @@ Improvements
73 73 New features
74 74 ------------
75 75  
  76 +* [`3606 <https://cta-redmine.irap.omp.eu/issues/3606>`_] -
  77 + Add ``GModelSpectralBins`` class
76 78 * [`3527 <https://cta-redmine.irap.omp.eu/issues/3527>`_] -
77 79 Add ``GCTABackground2D`` class
78 80 * [`3469 <https://cta-redmine.irap.omp.eu/issues/3469>`_] -
... ...
include/GModelSpectralBins.hpp 0 → 100644
  View file @8e14f9e
  1 +/***************************************************************************
  2 + * GModelSpectralBins.hpp - Spectral bins model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2021 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 GModelSpectralBins.hpp
  23 + * @brief Spectral bins model class definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GMODELSPECTRALBINS_HPP
  28 +#define GMODELSPECTRALBINS_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GModelSpectral.hpp"
  34 +#include "GModelPar.hpp"
  35 +#include "GEnergy.hpp"
  36 +
  37 +/* __ Forward declarations _______________________________________________ */
  38 +class GRan;
  39 +class GTime;
  40 +class GEbounds;
  41 +class GXmlElement;
  42 +
  43 +
  44 +/***********************************************************************//**
  45 + * @class GModelSpectralBins
  46 + *
  47 + * @brief Spectral bins model class
  48 + *
  49 + * This class implements spectral bins that have a constant intensity within
  50 + * their boundaries. A spectral bin is defined by a lower and upper energy
  51 + * limit, with the lower limit included and the upper limit excluded from an
  52 + * energy bin.
  53 + ***************************************************************************/
  54 +class GModelSpectralBins : public GModelSpectral {
  55 +
  56 +public:
  57 + // Constructors and destructors
  58 + GModelSpectralBins(void);
  59 + GModelSpectralBins(const GModelSpectral& model,
  60 + const GEbounds& ebounds,
  61 + const double& index = -2.0);
  62 + explicit GModelSpectralBins(const GXmlElement& xml);
  63 + GModelSpectralBins(const GModelSpectralBins& model);
  64 + virtual ~GModelSpectralBins(void);
  65 +
  66 + // Operators
  67 + virtual GModelSpectralBins& operator=(const GModelSpectralBins& model);
  68 +
  69 + // Implemented pure virtual base class methods
  70 + virtual void clear(void);
  71 + virtual GModelSpectralBins* clone(void) const;
  72 + virtual std::string classname(void) const;
  73 + virtual std::string type(void) const;
  74 + virtual double eval(const GEnergy& srcEng,
  75 + const GTime& srcTime = GTime(),
  76 + const bool& gradients = false) const;
  77 + virtual double flux(const GEnergy& emin,
  78 + const GEnergy& emax) const;
  79 + virtual double eflux(const GEnergy& emin,
  80 + const GEnergy& emax) const;
  81 + virtual GEnergy mc(const GEnergy& emin,
  82 + const GEnergy& emax,
  83 + const GTime& time,
  84 + GRan& ran) const;
  85 + virtual void read(const GXmlElement& xml);
  86 + virtual void write(GXmlElement& xml) const;
  87 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  88 +
  89 + // Other methods
  90 + int bins(void) const;
  91 + void append(const GEnergy& emin,
  92 + const GEnergy& emax,
  93 + const double& intensity);
  94 + void insert(const int& index,
  95 + const GEnergy& emin,
  96 + const GEnergy& emax,
  97 + const double& intensity);
  98 + void remove(const int& index);
  99 + void reserve(const int& num);
  100 + void extend(const GModelSpectralBins& bins);
  101 + double index(void) const;
  102 + void index(const double& index);
  103 + GEnergy emin(const int& index) const;
  104 + GEnergy emax(const int& index) const;
  105 + void emin(const int& index, const GEnergy& emin);
  106 + void emax(const int& index, const GEnergy& emax);
  107 + double intensity(const int& index) const;
  108 + void intensity(const int& index, const double& intensity);
  109 +
  110 +protected:
  111 + // Protected methods
  112 + void init_members(void);
  113 + void copy_members(const GModelSpectralBins& model);
  114 + void free_members(void);
  115 + int bin_index(const GEnergy& energy) const;
  116 + void update_pars(void);
  117 + void set_cache(void) const;
  118 + void mc_update(const GEnergy& emin, const GEnergy& emax) const;
  119 +
  120 + // Protected members
  121 + GModelPar m_index; //!< Spectral index of all bins
  122 + std::vector<GModelPar> m_emin; //!< Lower energy limits
  123 + std::vector<GModelPar> m_emax; //!< Upper energy limits
  124 + std::vector<GModelPar> m_values; //!< Bin values
  125 +
  126 + // Evaluation cache
  127 + mutable std::vector<double> m_epivot; //!< Power-law pivot energies in MeV
  128 +
  129 + // MC cache
  130 + mutable GEnergy m_mc_emin; //!< Minimum energy
  131 + mutable GEnergy m_mc_emax; //!< Maximum energy
  132 + mutable std::vector<double> m_mc_cum; //!< Cumulative distribution
  133 + mutable std::vector<double> m_mc_min; //!< Lower boundary for MC
  134 + mutable std::vector<double> m_mc_max; //!< Upper boundary for MC
  135 + mutable std::vector<double> m_mc_exp; //!< Exponent for MC
  136 +};
  137 +
  138 +
  139 +/***********************************************************************//**
  140 + * @brief Return class name
  141 + *
  142 + * @return String containing the class name ("GModelSpectralBins").
  143 + ***************************************************************************/
  144 +inline
  145 +std::string GModelSpectralBins::classname(void) const
  146 +{
  147 + return ("GModelSpectralBins");
  148 +}
  149 +
  150 +
  151 +/***********************************************************************//**
  152 + * @brief Return model type
  153 + *
  154 + * @return "BinFunction".
  155 + *
  156 + * Returns the type of the spectral bin function model.
  157 + ***************************************************************************/
  158 +inline
  159 +std::string GModelSpectralBins::type(void) const
  160 +{
  161 + return "BinFunction";
  162 +}
  163 +
  164 +
  165 +/***********************************************************************//**
  166 + * @brief Return number of bins
  167 + *
  168 + * @return Number of bins.
  169 + *
  170 + * Returns the number of bins in the bin function model.
  171 + ***************************************************************************/
  172 +inline
  173 +int GModelSpectralBins::bins(void) const
  174 +{
  175 + return (int)m_values.size();
  176 +}
  177 +
  178 +#endif /* GMODELSPECTRALBINS_HPP */
... ...
include/GModelSpectralNodes.hpp
  View file @8e14f9e
... ... @@ -133,7 +133,6 @@ protected:
133 133 void copy_members(const GModelSpectralNodes& model);
134 134 void free_members(void);
135 135 void update_pars(void);
136   - void load_nodes(const std::string& filename);
137 136 void set_cache(void) const;
138 137 void set_eval_cache(void) const;
139 138 void set_flux_cache(void) const;
... ...
include/GammaLib.hpp
  View file @8e14f9e
1 1 /***************************************************************************
2 2 * GammaLib.hpp - Gamma-Ray Astronomy Library Header file *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2006-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2006-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -241,6 +241,7 @@
241 241 #include "GModelSpectralGauss.hpp"
242 242 #include "GModelSpectralLogParabola.hpp"
243 243 #include "GModelSpectralMultiplicative.hpp"
  244 +#include "GModelSpectralBins.hpp"
244 245 #include "GModelSpectralNodes.hpp"
245 246 #include "GModelSpectralPlaw.hpp"
246 247 #include "GModelSpectralPlawPhotonFlux.hpp"
... ...
include/Makefile.am
  View file @8e14f9e
1 1 #############################################################################
2 2 # Makefile.am for include module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2020 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2021 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -200,6 +200,7 @@ pkginclude_HEADERS = GammaLib.hpp \
200 200 GModelSpectralGauss.hpp \
201 201 GModelSpectralLogParabola.hpp \
202 202 GModelSpectralMultiplicative.hpp \
  203 + GModelSpectralBins.hpp \
203 204 GModelSpectralNodes.hpp \
204 205 GModelSpectralPlaw.hpp \
205 206 GModelSpectralPlawPhotonFlux.hpp \
... ...
pyext/GModelSpectralBins.i 0 → 100644
  View file @8e14f9e
  1 +/***************************************************************************
  2 + * GModelSpectralBins.i - Spectral bins model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2021 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 GModelSpectralBins.i
  23 + * @brief Spectral bins model class Python interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelSpectralBins.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GModelSpectralBins
  34 + *
  35 + * @brief Spectral bins model class
  36 + ***************************************************************************/
  37 +class GModelSpectralBins : public GModelSpectral {
  38 +public:
  39 + // Constructors and destructors
  40 + GModelSpectralBins(void);
  41 + GModelSpectralBins(const GModelSpectral& model,
  42 + const GEbounds& ebounds,
  43 + const double& index = -2.0);
  44 + explicit GModelSpectralBins(const GXmlElement& xml);
  45 + GModelSpectralBins(const GModelSpectralBins& model);
  46 + virtual ~GModelSpectralBins(void);
  47 +
  48 + // Implemented pure virtual methods
  49 + virtual void clear(void);
  50 + virtual GModelSpectralBins* clone(void) const;
  51 + virtual std::string classname(void) const;
  52 + virtual std::string type(void) const;
  53 + virtual double eval(const GEnergy& srcEng,
  54 + const GTime& srcTime = GTime(),
  55 + const bool& gradients = false) const;
  56 + virtual double flux(const GEnergy& emin,
  57 + const GEnergy& emax) const;
  58 + virtual double eflux(const GEnergy& emin,
  59 + const GEnergy& emax) const;
  60 + virtual GEnergy mc(const GEnergy& emin,
  61 + const GEnergy& emax,
  62 + const GTime& time,
  63 + GRan& ran) const;
  64 + virtual void read(const GXmlElement& xml);
  65 + virtual void write(GXmlElement& xml) const;
  66 +
  67 + // Other methods
  68 + int bins(void) const;
  69 + void append(const GEnergy& emin,
  70 + const GEnergy& emax,
  71 + const double& intensity);
  72 + void insert(const int& index,
  73 + const GEnergy& emin,
  74 + const GEnergy& emax,
  75 + const double& intensity);
  76 + void remove(const int& index);
  77 + void reserve(const int& num);
  78 + void extend(const GModelSpectralBins& bins);
  79 + GEnergy emin(const int& index) const;
  80 + GEnergy emax(const int& index) const;
  81 + void emin(const int& index, const GEnergy& emin);
  82 + void emax(const int& index, const GEnergy& emax);
  83 + double intensity(const int& index) const;
  84 + void intensity(const int& index, const double& intensity);
  85 +};
  86 +
  87 +
  88 +/***********************************************************************//**
  89 + * @brief GModelSpectralBins class extension
  90 + ***************************************************************************/
  91 +%extend GModelSpectralBins {
  92 + GModelSpectralBins copy() {
  93 + return (*self);
  94 + }
  95 +%pythoncode {
  96 + def __getstate__(self):
  97 + xml = gammalib.GXmlElement()
  98 + self.write(xml)
  99 + state = xml,
  100 + return state
  101 + def __setstate__(self, state):
  102 + if state[0].elements('bin') == 0:
  103 + self.__init__()
  104 + else:
  105 + self.__init__(state[0])
  106 +}
  107 +};
... ...
pyext/gammalib/model.i
  View file @8e14f9e
1 1 /***************************************************************************
2 2 * model.i - Model module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2020 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2021 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -206,6 +206,7 @@
206 206 %include "GModelSpectralFunc.i"
207 207 %include "GModelSpectralGauss.i"
208 208 %include "GModelSpectralLogParabola.i"
  209 +%include "GModelSpectralBins.i"
209 210 %include "GModelSpectralNodes.i"
210 211 %include "GModelSpectralPlaw.i"
211 212 %include "GModelSpectralPlawPhotonFlux.i"
... ...
src/model/GModelSpectralBins.cpp 0 → 100644
  View file @8e14f9e
  1 +/***************************************************************************
  2 + * GModelSpectralBins.cpp - Spectral bins model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2021 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 GModelSpectralBins.cpp
  23 + * @brief Spectral bins model class implementation
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cmath>
  32 +#include "GException.hpp"
  33 +#include "GTools.hpp"
  34 +#include "GRan.hpp"
  35 +#include "GEbounds.hpp"
  36 +#include "GModelSpectralBins.hpp"
  37 +#include "GModelSpectralRegistry.hpp"
  38 +
  39 +/* __ Constants __________________________________________________________ */
  40 +
  41 +/* __ Globals ____________________________________________________________ */
  42 +const GModelSpectralBins g_spectral_bins_seed;
  43 +const GModelSpectralRegistry g_spectral_bins_registry(&g_spectral_bins_seed);
  44 +
  45 +/* __ Method name definitions ____________________________________________ */
  46 +#define G_FLUX "GModelSpectralBins::flux(GEnergy&, GEnergy&)"
  47 +#define G_EFLUX "GModelSpectralBins::eflux(GEnergy&, GEnergy&)"
  48 +#define G_MC "GModelSpectralBins::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
  49 +#define G_READ "GModelSpectralBins::read(GXmlElement&)"
  50 +#define G_WRITE "GModelSpectralBins::write(GXmlElement&)"
  51 +#define G_APPEND "GModelSpectralBins::append(GEnergy&, GEnergy&, double&)"
  52 +#define G_INSERT "GModelSpectralBins::insert(int&, GEnergy&, GEnergy&, "\
  53 + "double&)"
  54 +#define G_REMOVE "GModelSpectralBins::remove(int&)"
  55 +#define G_EMIN_GET "GModelSpectralBins::emin(int&)"
  56 +#define G_EMAX_GET "GModelSpectralBins::emax(int&)"
  57 +#define G_EMIN_SET "GModelSpectralBins::emin(int&, GEnergy&)"
  58 +#define G_EMAX_SET "GModelSpectralBins::emax(int&, GEnergy&)"
  59 +#define G_INTENSITY_GET "GModelSpectralBins::intensity(int&)"
  60 +#define G_INTENSITY_SET "GModelSpectralBins::intensity(int&, double&)"
  61 +
  62 +/* __ Macros _____________________________________________________________ */
  63 +
  64 +/* __ Coding definitions _________________________________________________ */
  65 +
  66 +/* __ Debug definitions __________________________________________________ */
  67 +
  68 +
  69 +/*==========================================================================
  70 + = =
  71 + = Constructors/destructors =
  72 + = =
  73 + ==========================================================================*/
  74 +
  75 +/***********************************************************************//**
  76 + * @brief Void constructor
  77 + *
  78 + * Constructs an empty spectral bin model.
  79 + ***************************************************************************/
  80 +GModelSpectralBins::GModelSpectralBins(void) : GModelSpectral()
  81 +{
  82 + // Initialise members
  83 + init_members();
  84 +
  85 + // Return
  86 + return;
  87 +}
  88 +
  89 +
  90 +/***********************************************************************//**
  91 + * @brief Spectral model constructor
  92 + *
  93 + * @param[in] model Spectral model.
  94 + * @param[in] ebounds Energy boundaries.
  95 + * @param[in] index Spectral index.
  96 + *
  97 + * Constructs a spectral bins model from any spectral model using the
  98 + * energy boundaries that are specified by @p ebounds. Within each spectral
  99 + * bin the intensity follows a power law with spectral index @p index.
  100 + ***************************************************************************/
  101 +GModelSpectralBins::GModelSpectralBins(const GModelSpectral& model,
  102 + const GEbounds& ebounds,
  103 + const double& index) :
  104 + GModelSpectral()
  105 +{
  106 + // Initialise members
  107 + init_members();
  108 +
  109 + // Append bins for all energies
  110 + for (int i = 0; i < ebounds.size(); ++i) {
  111 + append(ebounds.emin(i), ebounds.emax(i), model.eval(ebounds.elogmean(i)));
  112 + }
  113 +
  114 + // Store index
  115 + this->index(index);
  116 +
  117 + // Return
  118 + return;
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief XML constructor
  124 + *
  125 + * @param[in] xml XML element.
  126 + *
  127 + * Construct spectral bins model by extracting information from an XML
  128 + * element. See the read() method for more information about the expected
  129 + * structure of the XML element.
  130 + ***************************************************************************/
  131 +GModelSpectralBins::GModelSpectralBins(const GXmlElement& xml) :
  132 + GModelSpectral()
  133 +{
  134 + // Initialise members
  135 + init_members();
  136 +
  137 + // Read information from XML element
  138 + read(xml);
  139 +
  140 + // Return
  141 + return;
  142 +}
  143 +
  144 +
  145 +/***********************************************************************//**
  146 + * @brief Copy constructor
  147 + *
  148 + * @param[in] model Spectral bins model.
  149 + ***************************************************************************/
  150 +GModelSpectralBins::GModelSpectralBins(const GModelSpectralBins& model) :
  151 + GModelSpectral(model)
  152 +{
  153 + // Initialise members
  154 + init_members();
  155 +
  156 + // Copy members
  157 + copy_members(model);
  158 +
  159 + // Return
  160 + return;
  161 +}
  162 +
  163 +
  164 +/***********************************************************************//**
  165 + * @brief Destructor
  166 + ***************************************************************************/
  167 +GModelSpectralBins::~GModelSpectralBins(void)
  168 +{
  169 + // Free members
  170 + free_members();
  171 +
  172 + // Return
  173 + return;
  174 +}
  175 +
  176 +
  177 +/*==========================================================================
  178 + = =
  179 + = Operators =
  180 + = =
  181 + ==========================================================================*/
  182 +
  183 +/***********************************************************************//**
  184 + * @brief Assignment operator
  185 + *
  186 + * @param[in] model Spectral bins model.
  187 + * @return Spectral bins model.
  188 + ***************************************************************************/
  189 +GModelSpectralBins& GModelSpectralBins::operator=(const GModelSpectralBins& model)
  190 +{
  191 + // Execute only if object is not identical
  192 + if (this != &model) {
  193 +
  194 + // Copy base class members
  195 + this->GModelSpectral::operator=(model);
  196 +
  197 + // Free members
  198 + free_members();
  199 +
  200 + // Initialise members
  201 + init_members();
  202 +
  203 + // Copy members
  204 + copy_members(model);
  205 +
  206 + } // endif: object was not identical
  207 +
  208 + // Return
  209 + return *this;
  210 +}
  211 +
  212 +
  213 +/*==========================================================================
  214 + = =
  215 + = Public methods =
  216 + = =
  217 + ==========================================================================*/
  218 +
  219 +/***********************************************************************//**
  220 + * @brief Clear spectral bins model
  221 + ***************************************************************************/
  222 +void GModelSpectralBins::clear(void)
  223 +{
  224 + // Free class members (base and derived classes, derived class first)
  225 + free_members();
  226 + this->GModelSpectral::free_members();
  227 +
  228 + // Initialise members
  229 + this->GModelSpectral::init_members();
  230 + init_members();
  231 +
  232 + // Return
  233 + return;
  234 +}
  235 +
  236 +
  237 +/***********************************************************************//**
  238 + * @brief Clone spectral bins model
  239 +***************************************************************************/
  240 +GModelSpectralBins* GModelSpectralBins::clone(void) const
  241 +{
  242 + // Clone spectral bins model
  243 + return new GModelSpectralBins(*this);
  244 +}
  245 +
  246 +
  247 +/***********************************************************************//**
  248 + * @brief Evaluate function
  249 + *
  250 + * @param[in] srcEng True photon energy.
  251 + * @param[in] srcTime True photon arrival time.
  252 + * @param[in] gradients Compute gradients?
  253 + * @return Model value (ph/cm2/s/MeV).
  254 + *
  255 + * @todo Document method.
  256 + ***************************************************************************/
  257 +double GModelSpectralBins::eval(const GEnergy& srcEng,
  258 + const GTime& srcTime,
  259 + const bool& gradients) const
  260 +{
  261 + // Initialise value
  262 + double value = 0.0;
  263 +
  264 + // Optionally initialise partial derivatives
  265 + if (gradients) {
  266 + m_index.factor_gradient(0.0);
  267 + for (int i = 0; i < m_values.size(); ++i) {
  268 + m_values[i].factor_gradient(0.0);
  269 + }
  270 + }
  271 +
  272 + // Get bin index
  273 + int ibin = bin_index(srcEng);
  274 +
  275 + // Continue only if index is valid
  276 + if (ibin >= 0) {
  277 +
  278 + // Get power-law index
  279 + double index = m_index.value();
  280 +
  281 + // Compute and store value
  282 + double eng = srcEng.MeV();
  283 + double e_norm = eng / m_epivot[ibin];
  284 + double power = (index != 0.0) ? std::pow(e_norm, index) : 1.0;
  285 +
  286 + // Compute function value
  287 + value = m_values[ibin].value() * power;
  288 +
  289 + // Compute gradients
  290 + if (gradients) {
  291 +
  292 + // Compute partial derivatives of the parameter values
  293 + double g_norm = (m_values[ibin].is_free())
  294 + ? m_values[ibin].scale() * power : 0.0;
  295 + double g_index = (m_index.is_free())
  296 + ? value * m_index.scale() * std::log(e_norm) : 0.0;
  297 +
  298 + // Set gradients
  299 + m_values[ibin].factor_gradient(g_norm);
  300 + m_index.factor_gradient(g_index);
  301 +
  302 + } // endif: gradient computation was requested
  303 +
  304 + // Compile option: Check for NaN/Inf
  305 + #if defined(G_NAN_CHECK)
  306 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  307 + std::cout << "*** ERROR: GModelSpectralBins::eval";
  308 + std::cout << "(srcEng=" << srcEng;
  309 + std::cout << ", srcTime=" << srcTime << "):";
  310 + std::cout << " NaN/Inf encountered";
  311 + std::cout << " (value=" << value;
  312 + std::cout << ", power=" << power;
  313 + std::cout << ")" << std::endl;
  314 + }
  315 + #endif
  316 +
  317 + } // endif: bin index was valid
  318 +
  319 + // Return
  320 + return value;
  321 +}
  322 +
  323 +
  324 +/***********************************************************************//**
  325 + * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
  326 + *
  327 + * @param[in] emin Minimum photon energy.
  328 + * @param[in] emax Maximum photon energy.
  329 + * @return Photon flux (ph/cm2/s).
  330 + *
  331 + * Computes
  332 + *
  333 + * \f[
  334 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
  335 + * \f]
  336 + *
  337 + * where
  338 + * - [@p emin, @p emax] is an energy interval, and
  339 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  340 + ***************************************************************************/
  341 +double GModelSpectralBins::flux(const GEnergy& emin,
  342 + const GEnergy& emax) const
  343 +{
  344 + // Initialise flux
  345 + double flux = 0.0;
  346 +
  347 + // Compute only if integration range is valid
  348 + if (emin < emax) {
  349 +
  350 + // Get energy range in MeV
  351 + double e_min = emin.MeV();
  352 + double e_max = emax.MeV();
  353 +
  354 + // Loop over all bins
  355 + for (int ibin = 0; ibin < bins(); ++ibin) {
  356 +
  357 + // Get energy boundaires of bin in MeV
  358 + double e_bin_min = m_emin[ibin].value();
  359 + double e_bin_max = m_emax[ibin].value();
  360 +
  361 + // Skip bins that do not overlap with energy range
  362 + if (e_min >= e_bin_max || e_max < e_bin_min) {
  363 + continue;
  364 + }
  365 +
  366 + // Raise lower boundary to minimum energy
  367 + if (e_bin_min < e_min) {
  368 + e_bin_min = e_min;
  369 + }
  370 +
  371 + // Lower upper boundary to maximum energy
  372 + if (e_bin_max > e_max) {
  373 + e_bin_max = e_max;
  374 + }
  375 +
  376 + // If energy interval is positive then compute flux
  377 + if (e_bin_max > e_bin_min) {
  378 +
  379 + // Add photon flux
  380 + flux += m_values[ibin].value() *
  381 + gammalib::plaw_photon_flux(e_bin_min,
  382 + e_bin_max,
  383 + m_epivot[ibin],
  384 + m_index.value());
  385 +
  386 + } // endif: energy interval was positive
  387 +
  388 + } // endfor: looped over all spectral bins
  389 +
  390 + } // endif: integration range was valid
  391 +
  392 + // Return
  393 + return flux;
  394 +}
  395 +
  396 +
  397 +/***********************************************************************//**
  398 + * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
  399 + *
  400 + * @param[in] emin Minimum photon energy.
  401 + * @param[in] emax Maximum photon energy.
  402 + * @return Energy flux (erg/cm2/s).
  403 + *
  404 + * Computes
  405 + *
  406 + * \f[
  407 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
  408 + * \f]
  409 + *
  410 + * where
  411 + * - [@p emin, @p emax] is an energy interval, and
  412 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  413 + ***************************************************************************/
  414 +double GModelSpectralBins::eflux(const GEnergy& emin,
  415 + const GEnergy& emax) const
  416 +{
  417 + // Initialise flux
  418 + double eflux = 0.0;
  419 +
  420 + // Compute only if integration range is valid
  421 + if (emin < emax) {
  422 +
  423 +
  424 + // Get energy range in MeV
  425 + double e_min = emin.MeV();
  426 + double e_max = emax.MeV();
  427 +
  428 + // Loop over all bins
  429 + for (int ibin = 0; ibin < bins(); ++ibin) {
  430 +
  431 + // Get energy boundaires of bin in MeV
  432 + double e_bin_min = m_emin[ibin].value();
  433 + double e_bin_max = m_emax[ibin].value();
  434 +
  435 + // Skip bins that do not overlap with energy range
  436 + if (e_min >= e_bin_max || e_max < e_bin_min) {
  437 + continue;
  438 + }
  439 +
  440 + // Raise lower boundary to minimum energy
  441 + if (e_bin_min < e_min) {
  442 + e_bin_min = e_min;
  443 + }
  444 +
  445 + // Lower upper boundary to maximum energy
  446 + if (e_bin_max > e_max) {
  447 + e_bin_max = e_max;
  448 + }
  449 +
  450 + // If energy interval is positive then compute flux
  451 + if (e_bin_max > e_bin_min) {
  452 +
  453 + // Add energy flux
  454 + eflux += m_values[ibin].value() *
  455 + gammalib::plaw_energy_flux(e_bin_min,
  456 + e_bin_max,
  457 + m_epivot[ibin],
  458 + m_index.value());
  459 +
  460 + } // endif: energy interval was positive
  461 +
  462 + } // endfor: looped over all spectral bins
  463 +
  464 + } // endif: integration range was valid
  465 +
  466 + // Return flux
  467 + return eflux;
  468 +}
  469 +
  470 +
  471 +/***********************************************************************//**
  472 + * @brief Returns MC energy between [emin, emax]
  473 + *
  474 + * @param[in] emin Minimum photon energy.
  475 + * @param[in] emax Maximum photon energy.
  476 + * @param[in] time True photon arrival time.
  477 + * @param[in,out] ran Random number generator.
  478 + * @return Energy.
  479 + *
  480 + * @exception GException::erange_invalid
  481 + * Energy range is invalid (emin < emax required).
  482 + * @exception GException::invalid_return_value
  483 + * No valid Monte Carlo cache
  484 + *
  485 + * Returns Monte Carlo energy by randomly drawing from bin function.
  486 + ***************************************************************************/
  487 +GEnergy GModelSpectralBins::mc(const GEnergy& emin,
  488 + const GEnergy& emax,
  489 + const GTime& time,
  490 + GRan& ran) const
  491 +{
  492 + // Throw an exception if energy range is invalid
  493 + if (emin >= emax) {
  494 + throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
  495 + "Minimum energy < maximum energy required.");
  496 + }
  497 +
  498 + // Allocate energy
  499 + GEnergy energy;
  500 +
  501 + // Update cache
  502 + mc_update(emin, emax);
  503 +
  504 + // Determine in which bin we reside. If this fails then thrown an
  505 + // exception
  506 + int inx = 0;
  507 + if (m_mc_cum.size() > 1) {
  508 + double u = ran.uniform();
  509 + for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
  510 + if (m_mc_cum[inx-1] <= u) {
  511 + break;
  512 + }
  513 + }
  514 + }
  515 + else if (m_mc_cum.size() == 0) {
  516 + std::string msg = "No valid bins found for energy interval ["+
  517 + emin.print()+","+emax.print()+"]. Either restrict "
  518 + "the energy range to the one covered by the "
  519 + "spectral bins or extend the spectral bins "
  520 + "in energy.";
  521 + throw GException::invalid_return_value(G_MC, msg);
  522 + }
  523 +
  524 + // Get random energy for specific bin
  525 + if (m_mc_exp[inx] != 0.0) {
  526 + double e_min = m_mc_min[inx];
  527 + double e_max = m_mc_max[inx];
  528 + double u = ran.uniform();
  529 + double eng = (u > 0.0)
  530 + ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
  531 + : 0.0;
  532 + energy.MeV(eng);
  533 + }
  534 + else {
  535 + double e_min = m_mc_min[inx];
  536 + double e_max = m_mc_max[inx];
  537 + double u = ran.uniform();
  538 + double eng = std::exp(u * (e_max - e_min) + e_min);
  539 + energy.MeV(eng);
  540 + }
  541 +
  542 + // Return energy
  543 + return energy;
  544 +}
  545 +
  546 +
  547 +/***********************************************************************//**
  548 + * @brief Read model from XML element
  549 + *
  550 + * @param[in] xml XML element.
  551 + *
  552 + * @exception GException::invalid_value
  553 + * Lower energy limit larger than upper energy limit.
  554 + *
  555 + * Reads the spectral information from an XML element. The format of the XML
  556 + * elements is
  557 + *
  558 + * <spectrum type="BinFunction">
  559 + * <parameter name="Index" ../>
  560 + * <bin>
  561 + * <parameter name="LowerLimit" ../>
  562 + * <parameter name="UpperLimit" ../>
  563 + * <parameter name="Intensity" ../>
  564 + * </bin>
  565 + * ...
  566 + * <bin>
  567 + * <parameter name="LowerLimit" ../>
  568 + * <parameter name="UpperLimit" ../>
  569 + * <parameter name="Intensity" ../>
  570 + * </bin>
  571 + * </spectrum>
  572 + *
  573 + * @todo Check that bins are ordered
  574 + * @todo Check that energy boundaries are not overlapping
  575 + ***************************************************************************/
  576 +void GModelSpectralBins::read(const GXmlElement& xml)
  577 +{
  578 + // Free space for bins
  579 + m_emin.clear();
  580 + m_emax.clear();
  581 + m_values.clear();
  582 +
  583 + // Get index parameter pointer
  584 + const GXmlElement* index = gammalib::xml_get_par(G_READ, xml, m_index.name());
  585 +
  586 + // Read index parameter
  587 + m_index.read(*index);
  588 +
  589 + // Get number of bins from XML file
  590 + int bins = xml.elements("bin");
  591 +
  592 + // Loop over all bins
  593 + for (int i = 0; i < bins; ++i) {
  594 +
  595 + // Allocate bin parameters
  596 + GModelPar emin;
  597 + GModelPar emax;
  598 + GModelPar intensity;
  599 +
  600 + // Get bins
  601 + const GXmlElement* bin = xml.element("bin", i);
  602 +
  603 + // Get parameters
  604 + const GXmlElement* lpar = gammalib::xml_get_par(G_READ, *bin, "LowerLimit");
  605 + const GXmlElement* upar = gammalib::xml_get_par(G_READ, *bin, "UpperLimit");
  606 + const GXmlElement* ipar = gammalib::xml_get_par(G_READ, *bin, "Intensity");
  607 +
  608 + // Read parameters
  609 + emin.read(*lpar);
  610 + emax.read(*upar);
  611 + intensity.read(*ipar);
  612 +
  613 + // Throw an exception if either LowerLimit is larger than UpperLimit
  614 + if (emin.value() > emax.value()) {
  615 + std::string msg = "Lower energy limit "+gammalib::str(emin.value())+
  616 + " MeV is larger than upper energy limit "+
  617 + gammalib::str(emax.value())+" MeV. Please "
  618 + "specify lower energy limits that are smaller "
  619 + "than the upper energy limits.";
  620 + throw GException::invalid_value(G_READ, msg);
  621 + }
  622 +
  623 + // Push bin parameters on list
  624 + m_emin.push_back(emin);
  625 + m_emax.push_back(emax);
  626 + m_values.push_back(intensity);
  627 +
  628 + } // endfor: looped over bins
  629 +
  630 + // Update parameter mapping
  631 + update_pars();
  632 +
  633 + // Set pre-computation cache
  634 + set_cache();
  635 +
  636 + // Return
  637 + return;
  638 +}
  639 +
  640 +
  641 +/***********************************************************************//**
  642 + * @brief Write model into XML element
  643 + *
  644 + * @param[in] xml XML element into which model information is written.
  645 + *
  646 + * @exception GException::invalid_value
  647 + * Existing XML element is not of required type
  648 + * Invalid number of model parameters or bins found in XML element.
  649 + *
  650 + * Writes the spectral information into an XML element. The format of the XML
  651 + * element is
  652 + *
  653 + * <spectrum type="BinFunction">
  654 + * <parameter name="Index" ../>
  655 + * <bin>
  656 + * <parameter name="LowerLimit" ../>
  657 + * <parameter name="UpperLimit" ../>
  658 + * <parameter name="Intensity" ../>
  659 + * </bin>
  660 + * ...
  661 + * <bin>
  662 + * <parameter name="LowerLimit" ../>
  663 + * <parameter name="UpperLimit" ../>
  664 + * <parameter name="Intensity" ../>
  665 + * </bin>
  666 + * </spectrum>
  667 + ***************************************************************************/
  668 +void GModelSpectralBins::write(GXmlElement& xml) const
  669 +{
  670 + // Set model type
  671 + if (xml.attribute("type") == "") {
  672 + xml.attribute("type", type());
  673 + }
  674 +
  675 + // Verify model type
  676 + if (xml.attribute("type") != type()) {
  677 + std::string msg = "Spectral model type \""+xml.attribute("type")+"\" "
  678 + "is not of type \""+type()+"\".";
  679 + throw GException::invalid_value(G_WRITE, msg);
  680 + }
  681 +
  682 + // Get index parameter
  683 + GXmlElement* index = gammalib::xml_need_par(G_WRITE, xml, m_index.name());
  684 +
  685 + // Write index parameter
  686 + m_index.write(*index);
  687 +
  688 + // Determine number of bins
  689 + int bins = this->bins();
  690 +
  691 + // If XML element has 0 bins then append bins
  692 + if (xml.elements("bin") == 0) {
  693 + for (int i = 0; i < bins; ++i) {
  694 + xml.append(GXmlElement("bin"));
  695 + }
  696 + }
  697 +
  698 + // Verify that XML element has the required number of bins and elements.
  699 + // Recall that the power-law index also counts as one element.
  700 + if (xml.elements() != bins+1) {
  701 + std::string msg = "Spectral bins model contains "+
  702 + gammalib::str(xml.elements())+
  703 + " elements but "+gammalib::str(bins+1)+" elements "
  704 + "were expected.";
  705 + throw GException::invalid_value(G_WRITE, msg);
  706 + }
  707 + if (xml.elements("bin") != bins) {
  708 + std::string msg = "Spectral bins model contains "+
  709 + gammalib::str(xml.elements("bin"))+
  710 + " \"bin\" elements but "+gammalib::str(bins)+
  711 + " \"bin\" elements were expected.";
  712 + throw GException::invalid_value(G_WRITE, msg);
  713 + }
  714 +
  715 + // Loop over all bins
  716 + for (int i = 0; i < bins; ++i) {
  717 +
  718 + // Get bins
  719 + GXmlElement* bin = xml.element("bin", i);
  720 +
  721 + // Get copy of parameters so that we can change their names
  722 + GModelPar emin = m_emin[i];
  723 + GModelPar emax = m_emax[i];
  724 + GModelPar intensity = m_values[i];
  725 +
  726 + // Set names since we appended for internal usage the indices to the
  727 + // parameter names, and we want to get rid of them for writing the
  728 + // model into the XML element
  729 + emin.name("LowerLimit");
  730 + emax.name("UpperLimit");
  731 + intensity.name("Intensity");
  732 +
  733 + // Get XML parameters
  734 + GXmlElement* lpar = gammalib::xml_need_par(G_WRITE, *bin, emin.name());
  735 + GXmlElement* upar = gammalib::xml_need_par(G_WRITE, *bin, emax.name());
  736 + GXmlElement* ipar = gammalib::xml_need_par(G_WRITE, *bin, intensity.name());
  737 +
  738 + // Write parameters
  739 + emin.write(*lpar);
  740 + emax.write(*upar);
  741 + intensity.write(*ipar);
  742 +
  743 + } // endfor: looped over bins
  744 +
  745 + // Return
  746 + return;
  747 +}
  748 +
  749 +
  750 +/***********************************************************************//**
  751 + * @brief Append bin
  752 + *
  753 + * @param[in] emin Lower energy limit.
  754 + * @param[in] emax Upper energy limit.
  755 + * @param[in] intensity Bin intensity.
  756 + *
  757 + * @exception GException::invalid_argument
  758 + * Lower energy limit is larger than upper energy limit.
  759 + *
  760 + * Appends one bin to the bin function. By default the lower and upper energy
  761 + * limit of the bin is fixed while the intensity of the bin is free.
  762 + ***************************************************************************/
  763 +void GModelSpectralBins::append(const GEnergy& emin,
  764 + const GEnergy& emax,
  765 + const double& intensity)
  766 +{
  767 + // Throw an exception if lower limit is larger than upper limit
  768 + if (emin > emax) {
  769 + std::string msg = "Lower energy limit "+emin.print()+" is larger than "
  770 + "upper energy limit "+emax.print()+". Please "
  771 + "specify a lower energy limit that is smaller "
  772 + "than the upper energy limit.";
  773 + throw GException::invalid_value(G_APPEND, msg);
  774 + }
  775 +
  776 + // Allocate bin parameters
  777 + GModelPar par_emin;
  778 + GModelPar par_emax;
  779 + GModelPar par_intensity;
  780 +
  781 + // Set emin attributes
  782 + par_emin.name("LowerLimit");
  783 + par_emin.value(emin.MeV());
  784 + par_emin.unit("MeV");
  785 + par_emin.has_grad(false);
  786 + par_emin.fix();
  787 +
  788 + // Set emax attributes
  789 + par_emax.name("UpperLimit");
  790 + par_emax.value(emax.MeV());
  791 + par_emax.unit("MeV");
  792 + par_emax.has_grad(false);
  793 + par_emax.fix();
  794 +
  795 + // Set intensity attributes
  796 + par_intensity.name("Intensity");
  797 + par_intensity.value(intensity);
  798 + par_intensity.unit("ph/cm2/s/MeV");
  799 + par_intensity.has_grad(true);
  800 + par_intensity.free();
  801 +
  802 + // Append to bins
  803 + m_emin.push_back(par_emin);
  804 + m_emax.push_back(par_emax);
  805 + m_values.push_back(par_intensity);
  806 +
  807 + // Update parameter mapping
  808 + update_pars();
  809 +
  810 + // Set pre-computation cache
  811 + set_cache();
  812 +
  813 + // Return
  814 + return;
  815 +}
  816 +
  817 +
  818 +/***********************************************************************//**
  819 + * @brief Insert bin
  820 + *
  821 + * @param[in] index Bin index [0,...,bins()-1].
  822 + * @param[in] emin Lower energy limit.
  823 + * @param[in] emax Upper energy limit.
  824 + * @param[in] intensity Bin intensity.
  825 + *
  826 + * @exception GException::out_of_range
  827 + * Bin index is out of range.
  828 + * @exception GException::invalid_argument
  829 + * Lower energy limit is larger than upper energy limit.
  830 + *
  831 + * Inserts a bin into the bin function before the bin with the specified
  832 + * @p index. By default the lower and upper energy limit of the bin is fixed
  833 + * while the intensity of the bin is free.
  834 + ***************************************************************************/
  835 +void GModelSpectralBins::insert(const int& index,
  836 + const GEnergy& emin,
  837 + const GEnergy& emax,
  838 + const double& intensity)
  839 +{
  840 + // Raise exception if index is outside boundary
  841 + #if defined(G_RANGE_CHECK)
  842 + if (index < 0 || index >= bins()) {
  843 + throw GException::out_of_range(G_INSERT, "Bin index", index, bins());
  844 + }
  845 + #endif
  846 +
  847 + // Throw an exception if lower limit is larger than upper limit
  848 + if (emin > emax) {
  849 + std::string msg = "Lower energy limit "+emin.print()+" is larger than "
  850 + "upper energy limit "+emax.print()+". Please "
  851 + "specify a lower energy limit that is smaller "
  852 + "than the upper energy limit.";
  853 + throw GException::invalid_value(G_INSERT, msg);
  854 + }
  855 +
  856 + // Allocate bin parameters
  857 + GModelPar par_emin;
  858 + GModelPar par_emax;
  859 + GModelPar par_intensity;
  860 +
  861 + // Set emin attributes
  862 + par_emin.name("LowerLimit");
  863 + par_emin.value(emin.MeV());
  864 + par_emin.unit("MeV");
  865 + par_emin.has_grad(false);
  866 + par_emin.fix();
  867 +
  868 + // Set emax attributes
  869 + par_emax.name("UpperLimit");
  870 + par_emax.value(emax.MeV());
  871 + par_emax.unit("MeV");
  872 + par_emax.has_grad(false);
  873 + par_emax.fix();
  874 +
  875 + // Set intensity attributes
  876 + par_intensity.name("Intensity");
  877 + par_intensity.value(intensity);
  878 + par_intensity.unit("ph/cm2/s/MeV");
  879 + par_intensity.has_grad(true);
  880 + par_intensity.free();
  881 +
  882 + // Insert bin
  883 + m_emin.insert(m_emin.begin()+index, par_emin);
  884 + m_emax.insert(m_emax.begin()+index, par_emax);
  885 + m_values.insert(m_values.begin()+index, par_intensity);
  886 +
  887 + // Update parameter mapping
  888 + update_pars();
  889 +
  890 + // Set pre-computation cache
  891 + set_cache();
  892 +
  893 + // Return
  894 + return;
  895 +}
  896 +
  897 +
  898 +/***********************************************************************//**
  899 + * @brief Remove bin
  900 + *
  901 + * @param[in] index Bin index [0,...,bins()-1].
  902 + *
  903 + * @exception GException::out_of_range
  904 + * Bin index is out of range.
  905 + *
  906 + * Removes bin of specified @p index from the bin function.
  907 + ***************************************************************************/
  908 +void GModelSpectralBins::remove(const int& index)
  909 +{
  910 + // Raise exception if index is outside boundary
  911 + #if defined(G_RANGE_CHECK)
  912 + if (index < 0 || index >= bins()) {
  913 + throw GException::out_of_range(G_REMOVE, "Bin index", index, bins());
  914 + }
  915 + #endif
  916 +
  917 + // Erase energy and intensity
  918 + m_emin.erase(m_emin.begin() + index);
  919 + m_emax.erase(m_emax.begin() + index);
  920 + m_values.erase(m_values.begin() + index);
  921 +
  922 + // Update parameter mapping
  923 + update_pars();
  924 +
  925 + // Set pre-computation cache
  926 + set_cache();
  927 +
  928 + // Return
  929 + return;
  930 +}
  931 +
  932 +
  933 +/***********************************************************************//**
  934 + * @brief Reserve space for bins
  935 + *
  936 + * @param[in] num Number of reserved bins.
  937 + *
  938 + * Reserves space for @p num bins.
  939 + ***************************************************************************/
  940 +void GModelSpectralBins::reserve(const int& num)
  941 +{
  942 + // Reserve space
  943 + m_emin.reserve(num);
  944 + m_emax.reserve(num);
  945 + m_values.reserve(num);
  946 +
  947 + // Return
  948 + return;
  949 +}
  950 +
  951 +
  952 +/***********************************************************************//**
  953 + * @brief Append bins from bin function
  954 + *
  955 + * @param[in] bins Bin function.
  956 + *
  957 + * Appends all bins from a bin function to current object.
  958 + ***************************************************************************/
  959 +void GModelSpectralBins::extend(const GModelSpectralBins& bins)
  960 +{
  961 + // Get number of bins in bin function. Note that we extract the size
  962 + // first to avoid an endless loop that arises when a container is
  963 + // appended to itself.
  964 + int num = bins.bins();
  965 +
  966 + // Continue only if bin function is not empty
  967 + if (num > 0) {
  968 +
  969 + // Reserve enough space
  970 + reserve(this->bins() + num);
  971 +
  972 + // Loop over all bins and append them to the bin function
  973 + for (int i = 0; i < num; ++i) {
  974 + m_emin.push_back(bins.m_emin[i]);
  975 + m_emax.push_back(bins.m_emax[i]);
  976 + m_values.push_back(bins.m_values[i]);
  977 + }
  978 +
  979 + // Update parameter mapping
  980 + update_pars();
  981 +
  982 + // Set pre-computation cache
  983 + set_cache();
  984 +
  985 + } // endif: bin function was not empty
  986 +
  987 + // Return
  988 + return;
  989 +}
  990 +
  991 +
  992 +/***********************************************************************//**
  993 + * @brief Return spectral power law index
  994 + *
  995 + * @return Spectrel power law index.
  996 + *
  997 + * Returns the spectral power law index.
  998 + ***************************************************************************/
  999 +double GModelSpectralBins::index(void) const
  1000 +{
  1001 + // Return spectral power law index
  1002 + return (m_index.value());
  1003 +}
  1004 +
  1005 +
  1006 +/***********************************************************************//**
  1007 + * @brief Set spectral power law index
  1008 + *
  1009 + * @param[in] index Spectrel power law index.
  1010 + *
  1011 + * Set spectral power law index.
  1012 + ***************************************************************************/
  1013 +void GModelSpectralBins::index(const double& index)
  1014 +{
  1015 + // Set m_index
  1016 + m_index.value(index);
  1017 +
  1018 + // Return
  1019 + return;
  1020 +}
  1021 +
  1022 +
  1023 +/***********************************************************************//**
  1024 + * @brief Return lower energy limit of bin
  1025 + *
  1026 + * @param[in] index Bin index [0,...,bins()-1].
  1027 + * @return Lower energy limit of bin @p index.
  1028 + *
  1029 + * @exception GException::out_of_range
  1030 + * Index is out of range.
  1031 + *
  1032 + * Returns the lower energy limit of bin @p index.
  1033 + ***************************************************************************/
  1034 +GEnergy GModelSpectralBins::emin(const int& index) const
  1035 +{
  1036 + // Raise an exception if index is out of range
  1037 + #if defined(G_RANGE_CHECK)
  1038 + if (index < 0 || index >= bins()) {
  1039 + throw GException::out_of_range(G_EMIN_GET, "Bin index", index, bins());
  1040 + }
  1041 + #endif
  1042 +
  1043 + // Retrieve energy
  1044 + GEnergy energy;
  1045 + energy.MeV(m_emin[index].value());
  1046 +
  1047 + // Return energy
  1048 + return energy;
  1049 +}
  1050 +
  1051 +
  1052 +/***********************************************************************//**
  1053 + * @brief Return upper energy limit of bin
  1054 + *
  1055 + * @param[in] index Bin index [0,...,bins()-1].
  1056 + * @return Upper energy limit of bin @p index.
  1057 + *
  1058 + * @exception GException::out_of_range
  1059 + * Index is out of range.
  1060 + *
  1061 + * Returns the upper energy limit of bin @p index.
  1062 + ***************************************************************************/
  1063 +GEnergy GModelSpectralBins::emax(const int& index) const
  1064 +{
  1065 + // Raise an exception if index is out of range
  1066 + #if defined(G_RANGE_CHECK)
  1067 + if (index < 0 || index >= bins()) {
  1068 + throw GException::out_of_range(G_EMAX_GET, "Bin index", index, bins());
  1069 + }
  1070 + #endif
  1071 +
  1072 + // Retrieve energy
  1073 + GEnergy energy;
  1074 + energy.MeV(m_emax[index].value());
  1075 +
  1076 + // Return energy
  1077 + return energy;
  1078 +}
  1079 +
  1080 +
  1081 +/***********************************************************************//**
  1082 + * @brief Set lower energy limit of bin
  1083 + *
  1084 + * @param[in] index Bin index [0,...,bins()-1].
  1085 + * @param[in] emin Lower energy limit of bin @p index.
  1086 + *
  1087 + * @exception GException::out_of_range
  1088 + * Index is out of range.
  1089 + *
  1090 + * Sets the lower energy limit of bin @p index.
  1091 + ***************************************************************************/
  1092 +void GModelSpectralBins::emin(const int& index, const GEnergy& emin)
  1093 +{
  1094 + // Raise an exception if index is out of range
  1095 + #if defined(G_RANGE_CHECK)
  1096 + if (index < 0 || index >= bins()) {
  1097 + throw GException::out_of_range(G_EMIN_SET, "Bin index", index, bins());
  1098 + }
  1099 + #endif
  1100 +
  1101 + // Set energy
  1102 + m_emin[index].value(emin.MeV());
  1103 +
  1104 + // Set pre-computation cache
  1105 + set_cache();
  1106 +
  1107 + // Return
  1108 + return;
  1109 +}
  1110 +
  1111 +
  1112 +/***********************************************************************//**
  1113 + * @brief Set upper energy limit of bin
  1114 + *
  1115 + * @param[in] index Bin index [0,...,bins()-1].
  1116 + * @param[in] emax Upper energy limit of bin @p index.
  1117 + *
  1118 + * @exception GException::out_of_range
  1119 + * Index is out of range.
  1120 + *
  1121 + * Sets the upper energy limit of bin @p index.
  1122 + ***************************************************************************/
  1123 +void GModelSpectralBins::emax(const int& index, const GEnergy& emax)
  1124 +{
  1125 + // Raise an exception if index is out of range
  1126 + #if defined(G_RANGE_CHECK)
  1127 + if (index < 0 || index >= bins()) {
  1128 + throw GException::out_of_range(G_EMAX_SET, "Bin index", index, bins());
  1129 + }
  1130 + #endif
  1131 +
  1132 + // Set energy
  1133 + m_emax[index].value(emax.MeV());
  1134 +
  1135 + // Set pre-computation cache
  1136 + set_cache();
  1137 +
  1138 + // Return
  1139 + return;
  1140 +}
  1141 +
  1142 +
  1143 +/***********************************************************************//**
  1144 + * @brief Return bin intensity
  1145 + *
  1146 + * @param[in] index Bin index [0,...,bins()-1].
  1147 + * @return Intensity of bin @p index.
  1148 + *
  1149 + * @exception GException::out_of_range
  1150 + * Index is out of range.
  1151 + *
  1152 + * Returns the intensity of bin @p index.
  1153 + ***************************************************************************/
  1154 +double GModelSpectralBins::intensity(const int& index) const
  1155 +{
  1156 + // Raise an exception if index is out of range
  1157 + #if defined(G_RANGE_CHECK)
  1158 + if (index < 0 || index >= bins()) {
  1159 + throw GException::out_of_range(G_INTENSITY_GET, "Bin index", index, bins());
  1160 + }
  1161 + #endif
  1162 +
  1163 + // Return intensity
  1164 + return (m_values[index].value());
  1165 +}
  1166 +
  1167 +
  1168 +/***********************************************************************//**
  1169 + * @brief Set bin intensity
  1170 + *
  1171 + * @param[in] index Bin index [0,...,bins()-1].
  1172 + * @param[in] intensity Intensity of bin @p index.
  1173 + *
  1174 + * @exception GException::out_of_range
  1175 + * Index is out of range.
  1176 + *
  1177 + * Set the intensity of bin @p index.
  1178 + ***************************************************************************/
  1179 +void GModelSpectralBins::intensity(const int& index, const double& intensity)
  1180 +{
  1181 + // Raise an exception if index is out of range
  1182 + #if defined(G_RANGE_CHECK)
  1183 + if (index < 0 || index >= bins()) {
  1184 + throw GException::out_of_range(G_INTENSITY_SET, "Bin index", index, bins());
  1185 + }
  1186 + #endif
  1187 +
  1188 + // Set intensity
  1189 + m_values[index].value(intensity);
  1190 +
  1191 + // Return
  1192 + return;
  1193 +}
  1194 +
  1195 +
  1196 +/***********************************************************************//**
  1197 + * @brief Print bin function information
  1198 + *
  1199 + * @param[in] chatter Chattiness.
  1200 + * @return String containing bin function information.
  1201 + ***************************************************************************/
  1202 +std::string GModelSpectralBins::print(const GChatter& chatter) const
  1203 +{
  1204 + // Initialise result string
  1205 + std::string result;
  1206 +
  1207 + // Continue only if chatter is not silent
  1208 + if (chatter != SILENT) {
  1209 +
  1210 + // Append header
  1211 + result.append("=== GModelSpectralBins ===");
  1212 +
  1213 + // Append information
  1214 + result.append("\n"+gammalib::parformat("Number of bins"));
  1215 + result.append(gammalib::str(bins()));
  1216 + result.append("\n"+gammalib::parformat("Number of parameters"));
  1217 + result.append(gammalib::str(size()));
  1218 + for (int i = 0; i < size(); ++i) {
  1219 + result.append("\n"+m_pars[i]->print(chatter));
  1220 + }
  1221 +
  1222 + } // endif: chatter was not silent
  1223 +
  1224 + // Return result
  1225 + return result;
  1226 +}
  1227 +
  1228 +
  1229 +/*==========================================================================
  1230 + = =
  1231 + = Private methods =
  1232 + = =
  1233 + ==========================================================================*/
  1234 +
  1235 +/***********************************************************************//**
  1236 + * @brief Initialise class members
  1237 + ***************************************************************************/
  1238 +void GModelSpectralBins::init_members(void)
  1239 +{
  1240 + // Initialise powerlaw index
  1241 + m_index.clear();
  1242 + m_index.name("Index");
  1243 + m_index.scale(1.0);
  1244 + m_index.value(-2.0); // default: -2.0
  1245 + m_index.range(-10.0,+10.0); // range: [-10,+10]
  1246 + m_index.fix();
  1247 + m_index.gradient(0.0);
  1248 + m_index.has_grad(true);
  1249 +
  1250 + // Initialise bin parameter vectors
  1251 + m_emin.clear();
  1252 + m_emax.clear();
  1253 + m_values.clear();
  1254 +
  1255 + // Initialise evaluation cache
  1256 + m_epivot.clear();
  1257 +
  1258 + // Initialise MC cache
  1259 + m_mc_emin.clear();
  1260 + m_mc_emax.clear();
  1261 + m_mc_cum.clear();
  1262 + m_mc_min.clear();
  1263 + m_mc_max.clear();
  1264 + m_mc_exp.clear();
  1265 +
  1266 + // Update parameter mapping
  1267 + update_pars();
  1268 +
  1269 + // Return
  1270 + return;
  1271 +}
  1272 +
  1273 +
  1274 +/***********************************************************************//**
  1275 + * @brief Copy class members
  1276 + *
  1277 + * @param[in] model Spectral bin function.
  1278 + ***************************************************************************/
  1279 +void GModelSpectralBins::copy_members(const GModelSpectralBins& model)
  1280 +{
  1281 + // Copy bin energies and values
  1282 + m_index = model.m_index;
  1283 + m_emin = model.m_emin;
  1284 + m_emax = model.m_emax;
  1285 + m_values = model.m_values;
  1286 +
  1287 + // Copy evluation cache
  1288 + m_epivot = model.m_epivot;
  1289 +
  1290 + // Copy MC cache
  1291 + m_mc_emin = model.m_mc_emin;
  1292 + m_mc_emax = model.m_mc_emax;
  1293 + m_mc_cum = model.m_mc_cum;
  1294 + m_mc_min = model.m_mc_min;
  1295 + m_mc_max = model.m_mc_max;
  1296 + m_mc_exp = model.m_mc_exp;
  1297 +
  1298 + // Update parameter mapping
  1299 + update_pars();
  1300 +
  1301 + // Return
  1302 + return;
  1303 +}
  1304 +
  1305 +
  1306 +/***********************************************************************//**
  1307 + * @brief Delete class members
  1308 + ***************************************************************************/
  1309 +void GModelSpectralBins::free_members(void)
  1310 +{
  1311 + // Return
  1312 + return;
  1313 +}
  1314 +
  1315 +
  1316 +/***********************************************************************//**
  1317 + * @brief Return bin index for energy
  1318 + *
  1319 + * @param[in] energy Energy.
  1320 + * @return Bin index (-1 if energy bin was not found).
  1321 + ***************************************************************************/
  1322 +int GModelSpectralBins::bin_index(const GEnergy& energy) const
  1323 +{
  1324 + // Initialise bin index
  1325 + int bin = -1;
  1326 +
  1327 + // Get energy in MeV
  1328 + double energy_MeV = energy.MeV();
  1329 +
  1330 + // Search bin index
  1331 + for (int i = 0; i < bins(); ++i) {
  1332 + if (energy_MeV >= m_emin[i].value() && energy_MeV < m_emax[i].value()) {
  1333 + bin = i;
  1334 + break;
  1335 + }
  1336 + }
  1337 +
  1338 + // Return bin index
  1339 + return bin;
  1340 +}
  1341 +
  1342 +
  1343 +/***********************************************************************//**
  1344 + * @brief Update parameter mapping
  1345 + *
  1346 + * Sets the parameter pointers in the m_pars array, enabling iterating over
  1347 + * all model parameters. This method needs to be called after changing the
  1348 + * number of bins in the spectral model. The method needs not to be called
  1349 + * after value update.
  1350 + ***************************************************************************/
  1351 +void GModelSpectralBins::update_pars(void)
  1352 +{
  1353 + // Clear parameter pointer(s)
  1354 + m_pars.clear();
  1355 +
  1356 + // Push back common spectral index
  1357 + m_pars.push_back(&m_index);
  1358 +
  1359 + // Set parameter pointers for all bins
  1360 + for (int i = 0; i < bins(); ++i) {
  1361 +
  1362 + // Set parameter names
  1363 + std::string lowerlimit_name = "LowerLimit" + gammalib::str(i);
  1364 + std::string upperlimit_name = "UpperLimit" + gammalib::str(i);
  1365 + std::string intensity_name = "Intensity" + gammalib::str(i);
  1366 +
  1367 + // Set lower limit attributes
  1368 + m_emin[i].name(lowerlimit_name);
  1369 + m_emin[i].unit("MeV");
  1370 + m_emin[i].has_grad(false);
  1371 +
  1372 + // Set upper limit attributes
  1373 + m_emax[i].name(upperlimit_name);
  1374 + m_emax[i].unit("MeV");
  1375 + m_emax[i].has_grad(false);
  1376 +
  1377 + // Set intensity attributes
  1378 + m_values[i].name(intensity_name);
  1379 + m_values[i].unit("ph/cm2/s/MeV");
  1380 + m_values[i].has_grad(true);
  1381 +
  1382 + // Set pointer
  1383 + m_pars.push_back(&(m_emin[i]));
  1384 + m_pars.push_back(&(m_emax[i]));
  1385 + m_pars.push_back(&(m_values[i]));
  1386 +
  1387 + } // endfor: looped over bins
  1388 +
  1389 + // Return
  1390 + return;
  1391 +}
  1392 +
  1393 +
  1394 +/***********************************************************************//**
  1395 + * @brief Set cache
  1396 + *
  1397 + * Sets the computation cache. The method computes for the time being the
  1398 + * pivot energies for all bins in MeV.
  1399 + ***************************************************************************/
  1400 +void GModelSpectralBins::set_cache(void) const
  1401 +{
  1402 + // Clear cache
  1403 + m_epivot.clear();
  1404 +
  1405 + // Compute pivot energies for all bins
  1406 + for (int i = 0; i < bins(); ++i) {
  1407 + double epivot = std::sqrt(m_emin[i].value() * m_emax[i].value());
  1408 + m_epivot.push_back(epivot);
  1409 + }
  1410 +
  1411 + // Return
  1412 + return;
  1413 +}
  1414 +
  1415 +
  1416 +/***********************************************************************//**
  1417 + * @brief Set MC pre-computation cache
  1418 + *
  1419 + * @param[in] emin Minimum energy.
  1420 + * @param[in] emax Maximum energy.
  1421 + *
  1422 + * This method sets up an array of indices and the cumulative distribution
  1423 + * function needed for MC simulations.
  1424 + ***************************************************************************/
  1425 +void GModelSpectralBins::mc_update(const GEnergy& emin,
  1426 + const GEnergy& emax) const
  1427 +{
  1428 + // Check if we need to update the cache
  1429 + if (emin != m_mc_emin || emax != m_mc_emax) {
  1430 +
  1431 + // Store new energy interval
  1432 + m_mc_emin = emin;
  1433 + m_mc_emax = emax;
  1434 +
  1435 + // Initialise cache
  1436 + m_mc_cum.clear();
  1437 + m_mc_min.clear();
  1438 + m_mc_max.clear();
  1439 + m_mc_exp.clear();
  1440 +
  1441 + // Get energy range in MeV
  1442 + double e_min = emin.MeV();
  1443 + double e_max = emax.MeV();
  1444 +
  1445 + // Continue only if e_max > e_min
  1446 + if (e_max > e_min) {
  1447 +
  1448 + // Loop over all bins
  1449 + for (int ibin = 0; ibin < bins(); ++ibin) {
  1450 +
  1451 + // Get energy boundaires of bin in MeV
  1452 + double e_bin_min = m_emin[ibin].value();
  1453 + double e_bin_max = m_emax[ibin].value();
  1454 +
  1455 + // Skip bins that do not overlap with energy range
  1456 + if (e_min >= e_bin_max || e_max < e_bin_min) {
  1457 + continue;
  1458 + }
  1459 +
  1460 + // Raise lower boundary to minimum energy
  1461 + if (e_bin_min < e_min) {
  1462 + e_bin_min = e_min;
  1463 + }
  1464 +
  1465 + // Lower upper boundary to maximum energy
  1466 + if (e_bin_max > e_max) {
  1467 + e_bin_max = e_max;
  1468 + }
  1469 +
  1470 + // If energy interval is positive then compute flux
  1471 + if (e_bin_max > e_bin_min) {
  1472 +
  1473 + // Get index
  1474 + double index = m_index.value();
  1475 +
  1476 + // Compute flux
  1477 + double flux = m_values[ibin].value() *
  1478 + gammalib::plaw_photon_flux(e_bin_min,
  1479 + e_bin_max,
  1480 + m_epivot[ibin],
  1481 + index);
  1482 +
  1483 + // Push results in cache
  1484 + m_mc_cum.push_back(flux);
  1485 + m_mc_min.push_back(e_bin_min);
  1486 + m_mc_max.push_back(e_bin_max);
  1487 + m_mc_exp.push_back(index);
  1488 +
  1489 + } // endif: energy interval was positive
  1490 +
  1491 + } // endfor: looped over all spectral bins
  1492 +
  1493 + // Build cumulative distribution
  1494 + for (int i = 1; i < m_mc_cum.size(); ++i) {
  1495 + m_mc_cum[i] += m_mc_cum[i-1];
  1496 + }
  1497 + double norm = m_mc_cum[m_mc_cum.size()-1];
  1498 + if (norm > 0.0) {
  1499 + for (int i = 0; i < m_mc_cum.size(); ++i) {
  1500 + m_mc_cum[i] /= norm;
  1501 + }
  1502 + }
  1503 +
  1504 + // Set MC values
  1505 + for (int i = 0; i < m_mc_cum.size(); ++i) {
  1506 +
  1507 + // Compute exponent
  1508 + double exponent = m_mc_exp[i] + 1.0;
  1509 +
  1510 + // Exponent dependend computation
  1511 + if (std::abs(exponent) > 1.0e-11) {
  1512 +
  1513 + // If the exponent is too small then use lower energy
  1514 + // boundary
  1515 + if (exponent < -50.0) {
  1516 + m_mc_exp[i] = 0.0;
  1517 + m_mc_min[i] = std::log(m_mc_min[i]);
  1518 + m_mc_max[i] = m_mc_min[i];
  1519 + }
  1520 +
  1521 + // ... otherwise if exponent is too large then use
  1522 + // upper energy boundary
  1523 + else if (exponent > +50.0) {
  1524 + m_mc_exp[i] = 0.0;
  1525 + m_mc_min[i] = std::log(m_mc_max[i]);
  1526 + m_mc_max[i] = m_mc_min[i];
  1527 + }
  1528 +
  1529 + // ... otherwise use transformation formula
  1530 + else {
  1531 + m_mc_exp[i] = exponent;
  1532 + m_mc_min[i] = std::pow(m_mc_min[i], exponent);
  1533 + m_mc_max[i] = std::pow(m_mc_max[i], exponent);
  1534 + }
  1535 + }
  1536 + else {
  1537 + m_mc_exp[i] = 0.0;
  1538 + m_mc_min[i] = std::log(m_mc_min[i]);
  1539 + m_mc_max[i] = std::log(m_mc_max[i]);
  1540 + }
  1541 +
  1542 + } // endfor: set MC values
  1543 +
  1544 + } // endif: e_max > e_min
  1545 +
  1546 + } // endif: Update was required
  1547 +
  1548 + // Return
  1549 + return;
  1550 +}
... ...
src/model/Makefile.am
  View file @8e14f9e
1 1 #############################################################################
2 2 # Makefile.am for model module #
3 3 # ------------------------------------------------------------------------- #
4   -# Copyright (C) 2016-2020 Juergen Knoedlseder #
  4 +# Copyright (C) 2016-2021 Juergen Knoedlseder #
5 5 # ------------------------------------------------------------------------- #
6 6 # #
7 7 # This program is free software: you can redistribute it and/or modify #
... ... @@ -67,6 +67,7 @@ sources = GModelPar.cpp \
67 67 GModelSpectralGauss.cpp \
68 68 GModelSpectralLogParabola.cpp \
69 69 GModelSpectralMultiplicative.cpp \
  70 + GModelSpectralBins.cpp \
70 71 GModelSpectralNodes.cpp \
71 72 GModelSpectralPlaw.cpp \
72 73 GModelSpectralPlawPhotonFlux.cpp \
... ...
test/data/model_point_bins.xml 0 → 100644
  View file @8e14f9e
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Crab" type="PointSource">
  4 + <spectrum type="BinFunction">
  5 + <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="0"/>
  6 + <bin>
  7 + <parameter scale="1.0" name="LowerLimit" min="0.1" max="1.0e20" value="0.75" free="0"/>
  8 + <parameter scale="1.0" name="UpperLimit" min="0.1" max="1.0e20" value="1.0" free="0"/>
  9 + <parameter scale="1e-07" name="Intensity" min="1e-07" max="1000.0" value="1.0" free="1"/>
  10 + </bin>
  11 + <bin>
  12 + <parameter scale="1.0" name="LowerLimit" min="0.1" max="1.0e20" value="1.0" free="0"/>
  13 + <parameter scale="1.0" name="UpperLimit" min="0.1" max="1.0e20" value="3.0" free="0"/>
  14 + <parameter scale="1e-07" name="Intensity" min="1e-07" max="1000.0" value="0.5" free="1"/>
  15 + </bin>
  16 + <bin>
  17 + <parameter scale="1.0" name="LowerLimit" min="0.1" max="1.0e20" value="3.0" free="0"/>
  18 + <parameter scale="1.0" name="UpperLimit" min="0.1" max="1.0e20" value="10.0" free="0"/>
  19 + <parameter scale="1e-07" name="Intensity" min="1e-07" max="1000.0" value="0.1" free="1"/>
  20 + </bin>
  21 + <bin>
  22 + <parameter scale="1.0" name="LowerLimit" min="0.1" max="1.0e20" value="10.0" free="0"/>
  23 + <parameter scale="1.0" name="UpperLimit" min="0.1" max="1.0e20" value="30.0" free="0"/>
  24 + <parameter scale="1e-07" name="Intensity" min="1e-07" max="1000.0" value="0.05" free="1"/>
  25 + </bin>
  26 + </spectrum>
  27 + <spatialModel type="PointSource">
  28 + <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
  29 + <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
  30 + </spatialModel>
  31 + </source>
  32 +</source_library>
... ...
test/test_GModel.cpp
  View file @8e14f9e
... ... @@ -95,6 +95,7 @@ void TestGModel::set(void)
95 95 m_xml_model_point_smoothbplaw = datadir + "/model_point_smoothbplaw.xml";
96 96 m_xml_model_point_supeplaw = datadir + "/model_point_supeplaw.xml";
97 97 m_xml_model_point_logparabola = datadir + "/model_point_logparabola.xml";
  98 + m_xml_model_point_bins = datadir + "/model_point_bins.xml";
98 99 m_xml_model_point_nodes = datadir + "/model_point_nodes.xml";
99 100 m_xml_model_point_filefct = datadir + "/model_point_filefct.xml";
100 101 m_xml_model_point_table = datadir + "/model_point_table.xml";
... ... @@ -201,6 +202,8 @@ void TestGModel::set(void)
201 202 "Test GModelSpectralExponential");
202 203 append(static_cast<pfunction>(&TestGModel::test_nodes),
203 204 "Test GModelSpectralNodes");
  205 + append(static_cast<pfunction>(&TestGModel::test_bins),
  206 + "Test GModelSpectralBins");
204 207 append(static_cast<pfunction>(&TestGModel::test_filefct),
205 208 "Test GModelSpectralFunc");
206 209 append(static_cast<pfunction>(&TestGModel::test_table),
... ... @@ -2437,6 +2440,164 @@ void TestGModel::test_exponential(void)
2437 2440  
2438 2441  
2439 2442 /***********************************************************************//**
  2443 + * @brief Test GModelSpectralBins class
  2444 + ***************************************************************************/
  2445 +void TestGModel::test_bins(void)
  2446 +{
  2447 + // Test void constructor
  2448 + GModelSpectralBins model1;
  2449 + test_value(model1.type(), "BinFunction", "Check void model type");
  2450 + test_value(model1.bins(), 0, "Check number of bins in void model");
  2451 + test_value(model1.size(), 1, "Check number of parameters in void model");
  2452 +
  2453 + // Test model constructor
  2454 + GModelSpectralPlaw plaw(2.0, -2.1, GEnergy(100.0, "MeV"));
  2455 + GEbounds ebounds(3, GEnergy(1.0, "MeV"), GEnergy(1000.0, "MeV"));
  2456 + GModelSpectralBins model2(plaw, ebounds, -3.0);
  2457 + test_value(model2.bins(), 3, "Check number of bins in power-law constructed model");
  2458 + test_value(model2.size(), 10, "Check number of parameters in power-law constructed model");
  2459 + test_value(model2.index(), -3.0, "Check index in power-law constructed model");
  2460 + test_value(model2.emin(0).MeV(), 1.0, "Check emin[0] in power-law constructed model");
  2461 + test_value(model2.emax(0).MeV(), 10.0, "Check emax[0] in power-law constructed model");
  2462 + test_value(model2.intensity(0), 2825.075089, "Check intensity[0] in power-law constructed model");
  2463 + test_value(model2.emin(1).MeV(), 10.0, "Check emin[1] in power-law constructed model");
  2464 + test_value(model2.emax(1).MeV(), 100.0, "Check emax[1] in power-law constructed model");
  2465 + test_value(model2.intensity(1), 22.440369, "Check intensity[1] in power-law constructed model");
  2466 + test_value(model2.emin(2).MeV(), 100.0, "Check emin[2] in power-law constructed model");
  2467 + test_value(model2.emax(2).MeV(), 1000.0, "Check emax[2] in power-law constructed model");
  2468 + test_value(model2.intensity(2), 0.1782502, "Check intensity[2] in power-law constructed model");
  2469 + test_value(model2.eval(GEnergy(1.0, "MeV")), 89336.718430, "Check eval(1 MeV) in power-law constructed model");
  2470 + test_value(model2.eval(GEnergy(9.0, "MeV")), 122.546939, "Check eval(9 MeV) in power-law constructed model");
  2471 + test_value(model2.eval(GEnergy(76.0, "MeV")), 1.61655029, "Check eval(76 MeV) in power-law constructed model");
  2472 +
  2473 + // Test XML constructor
  2474 + GXml xml(m_xml_model_point_bins);
  2475 + GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
  2476 + GModelSpectralBins model3(*element);
  2477 + test_value(model3.bins(), 4, "Check number of bins in XML constructed model");
  2478 + test_value(model3.size(), 13, "Check number of parameters in XML constructed model");
  2479 + test_value(model3.index(), -2.48, "Check index in XML constructed model");
  2480 + test_value(model3.emin(0).MeV(), 0.75, "Check emin[0] in XML constructed model");
  2481 + test_value(model3.emax(0).MeV(), 1.0, "Check emax[0] in XML constructed model");
  2482 + test_value(model3.intensity(0), 1.0e-7, "Check intensity[0] in XML constructed model");
  2483 + test_value(model3.emin(1).MeV(), 1.0, "Check emin[1] in XML constructed model");
  2484 + test_value(model3.emax(1).MeV(), 3.0, "Check emax[1] in XML constructed model");
  2485 + test_value(model3.intensity(1), 0.5e-7, "Check intensity[1] in XML constructed model");
  2486 + test_value(model3.emin(2).MeV(), 3.0, "Check emin[2] in XML constructed model");
  2487 + test_value(model3.emax(2).MeV(), 10.0, "Check emax[2] in XML constructed model");
  2488 + test_value(model3.intensity(2), 0.1e-7, "Check intensity[2] in XML constructed model");
  2489 + test_value(model3.emin(3).MeV(), 10.0, "Check emin[3] in XML constructed model");
  2490 + test_value(model3.emax(3).MeV(), 30.0, "Check emax[3] in XML constructed model");
  2491 + test_value(model3.intensity(3), 0.05e-7, "Check intensity[3] in XML constructed model");
  2492 +
  2493 + // Test copy constructor
  2494 + GModelSpectralBins model4(model2);
  2495 + test_value(model4.bins(), 3, "Check number of bins in power-law constructed model");
  2496 + test_value(model4.size(), 10, "Check number of parameters in power-law constructed model");
  2497 + test_value(model4.index(), -3.0, "Check index in power-law constructed model");
  2498 + test_value(model4.emin(0).MeV(), 1.0, "Check emin[0] in power-law constructed model");
  2499 + test_value(model4.emax(0).MeV(), 10.0, "Check emax[0] in power-law constructed model");
  2500 + test_value(model4.intensity(0), 2825.075089, "Check intensity[0] in power-law constructed model");
  2501 + test_value(model4.emin(1).MeV(), 10.0, "Check emin[1] in power-law constructed model");
  2502 + test_value(model4.emax(1).MeV(), 100.0, "Check emax[1] in power-law constructed model");
  2503 + test_value(model4.intensity(1), 22.440369, "Check intensity[1] in power-law constructed model");
  2504 + test_value(model4.emin(2).MeV(), 100.0, "Check emin[2] in power-law constructed model");
  2505 + test_value(model4.emax(2).MeV(), 1000.0, "Check emax[2] in power-law constructed model");
  2506 + test_value(model4.intensity(2), 0.1782502, "Check intensity[2] in power-law constructed model");
  2507 +
  2508 + // Test XML write and read methods
  2509 + GXmlElement xml_element;
  2510 + model4.write(xml_element);
  2511 + GModelSpectralBins model5;
  2512 + model5.read(xml_element);
  2513 + test_value(model5.bins(), 3, "Check number of bins in power-law constructed model");
  2514 +
  2515 + // Test bin function manipulation
  2516 + GModelSpectralBins model6;
  2517 + model6.reserve(3);
  2518 + test_value(model6.bins(), 0, "Check number of bins after first manipulation");
  2519 + test_value(model6.size(), 1, "Check number of parameters after first manipulation");
  2520 + model6.append(GEnergy(1.0, "MeV"), GEnergy(10.0, "MeV"), 1.0);
  2521 + test_value(model6.bins(), 1, "Check number of bins after second manipulation");
  2522 + test_value(model6.size(), 4, "Check number of parameters after second manipulation");
  2523 + test_value(model6.emin(0).MeV(), 1.0, "Check emin[0] after second manipulation");
  2524 + test_value(model6.emax(0).MeV(), 10.0, "Check emax[0] after second manipulation");
  2525 + test_value(model6.intensity(0), 1.0, "Check intensity[0] after second manipulation");
  2526 + model6.append(GEnergy(10.0, "MeV"), GEnergy(100.0, "MeV"), 0.1);
  2527 + test_value(model6.bins(), 2, "Check number of bins after third manipulation");
  2528 + test_value(model6.size(), 7, "Check number of parameters after third manipulation");
  2529 + test_value(model6.emin(0).MeV(), 1.0, "Check emin[0] after third manipulation");
  2530 + test_value(model6.emax(0).MeV(), 10.0, "Check emax[0] after third manipulation");
  2531 + test_value(model6.intensity(0), 1.0, "Check intensity[0] after third manipulation");
  2532 + test_value(model6.emin(1).MeV(), 10.0, "Check emin[1] after third manipulation");
  2533 + test_value(model6.emax(1).MeV(), 100.0, "Check emax[1] after third manipulation");
  2534 + test_value(model6.intensity(1), 0.1, "Check intensity[1] after third manipulation");
  2535 + model6.remove(0);
  2536 + test_value(model6.bins(), 1, "Check number of bins after forth manipulation");
  2537 + test_value(model6.size(), 4, "Check number of parameters after forth manipulation");
  2538 + test_value(model6.emin(0).MeV(), 10.0, "Check emin[0] after forth manipulation");
  2539 + test_value(model6.emax(0).MeV(), 100.0, "Check emax[0] after forth manipulation");
  2540 + test_value(model6.intensity(0), 0.1, "Check intensity[0] after forth manipulation");
  2541 + model6.insert(0, GEnergy(0.9, "MeV"), GEnergy(9.9, "MeV"), 0.9);
  2542 + test_value(model6.bins(), 2, "Check number of bins after fifth manipulation");
  2543 + test_value(model6.size(), 7, "Check number of parameters after fifth manipulation");
  2544 + test_value(model6.emin(0).MeV(), 0.9, "Check emin[0] after fifth manipulation");
  2545 + test_value(model6.emax(0).MeV(), 9.9, "Check emax[0] after fifth manipulation");
  2546 + test_value(model6.intensity(0), 0.9, "Check intensity[0] after fifth manipulation");
  2547 + test_value(model6.emin(1).MeV(), 10.0, "Check emin[1] after fifth manipulation");
  2548 + test_value(model6.emax(1).MeV(), 100.0, "Check emax[1] after fifth manipulation");
  2549 + test_value(model6.intensity(1), 0.1, "Check intensity[1] after fifth manipulation");
  2550 + model6.extend(model6);
  2551 + test_value(model6.bins(), 4, "Check number of bins after sixth manipulation");
  2552 + test_value(model6.size(), 13, "Check number of parameters after sixth manipulation");
  2553 + test_value(model6.emin(0).MeV(), 0.9, "Check emin[0] after sixth manipulation");
  2554 + test_value(model6.emax(0).MeV(), 9.9, "Check emax[0] after sixth manipulation");
  2555 + test_value(model6.intensity(0), 0.9, "Check intensity[0] after sixth manipulation");
  2556 + test_value(model6.emin(1).MeV(), 10.0, "Check emin[1] after sixth manipulation");
  2557 + test_value(model6.emax(1).MeV(), 100.0, "Check emax[1] after sixth manipulation");
  2558 + test_value(model6.intensity(1), 0.1, "Check intensity[1] after sixth manipulation");
  2559 + test_value(model6.emin(2).MeV(), 0.9, "Check emin[2] after sixth manipulation");
  2560 + test_value(model6.emax(2).MeV(), 9.9, "Check emax[2] after sixth manipulation");
  2561 + test_value(model6.intensity(2), 0.9, "Check intensity[2] after sixth manipulation");
  2562 + test_value(model6.emin(3).MeV(), 10.0, "Check emin[3] after sixth manipulation");
  2563 + test_value(model6.emax(3).MeV(), 100.0, "Check emax[3] after sixth manipulation");
  2564 + test_value(model6.intensity(3), 0.1, "Check intensity[3] after sixth manipulation");
  2565 + model6.clear();
  2566 + test_value(model6.bins(), 0, "Check number of bins after seventh manipulation");
  2567 + test_value(model6.size(), 1, "Check number of parameters after seventh manipulation");
  2568 +
  2569 + // Test emin method
  2570 + model5.emin(0, GEnergy(0.1, "GeV"));
  2571 + test_value(model5.emin(0).GeV(), 0.1, "Check emin() method");
  2572 +
  2573 + // Test emax method
  2574 + model5.emax(0, GEnergy(0.2, "GeV"));
  2575 + test_value(model5.emax(0).GeV(), 0.2, "Check emax() method");
  2576 +
  2577 + // Test intensity method
  2578 + model5.intensity(0, 2.0e-7);
  2579 + test_value(model5.intensity(0), 2.0e-7, "Check intensity() method");
  2580 +
  2581 + // Test operator access
  2582 + const char* strarray[] = {"Index", "LowerLimit0", "UpperLimit0", "Intensity0",
  2583 + "LowerLimit0", "UpperLimit0", "Intensity1"};
  2584 + for (int i = 0; i < 7; ++i) {
  2585 + std::string keyname(strarray[i]);
  2586 + model5[keyname].remove_range(); // To allow setting of any value
  2587 + model5[keyname].value(2.1);
  2588 + model5[keyname].error(1.9);
  2589 + model5[keyname].gradient(0.8);
  2590 + test_value(model5[keyname].value(), 2.1);
  2591 + test_value(model5[keyname].error(), 1.9);
  2592 + test_value(model5[keyname].gradient(), 0.8);
  2593 + }
  2594 +
  2595 + // Exit test
  2596 + return;
  2597 +}
  2598 +
  2599 +
  2600 +/***********************************************************************//**
2440 2601 * @brief Test GModelSpectralNodes class
2441 2602 ***************************************************************************/
2442 2603 void TestGModel::test_nodes(void)
... ... @@ -2445,7 +2606,7 @@ void TestGModel::test_nodes(void)
2445 2606 GModelSpectralNodes model1;
2446 2607 test_value(model1.type(), "NodeFunction", "Check void model type");
2447 2608  
2448   - // Test node function manimulation
  2609 + // Test node function manipulation
2449 2610 GModelSpectralNodes model2;
2450 2611 model2.reserve(3);
2451 2612 test_value(model2.size(), 0);
... ...
test/test_GModel.hpp
  View file @8e14f9e
1 1 /***************************************************************************
2 2 * test_GModel.hpp - Test model module *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2020 by Jean-Baptiste Cayrou *
  4 + * copyright (C) 2012-2021 by Jean-Baptiste Cayrou *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -76,6 +76,7 @@ public:
76 76 void test_supeplaw(void);
77 77 void test_logparabola(void);
78 78 void test_multiplicative(void);
  79 + void test_bins(void);
79 80 void test_nodes(void);
80 81 void test_filefct(void);
81 82 void test_spectral_composite(void);
... ... @@ -130,6 +131,7 @@ private:
130 131 std::string m_xml_model_point_smoothbplaw;
131 132 std::string m_xml_model_point_supeplaw;
132 133 std::string m_xml_model_point_logparabola;
  134 + std::string m_xml_model_point_bins;
133 135 std::string m_xml_model_point_nodes;
134 136 std::string m_xml_model_point_filefct;
135 137 std::string m_xml_model_point_table;
... ...
test/test_GModel.py
  View file @8e14f9e
... ... @@ -143,6 +143,7 @@ class Test(gammalib.GPythonTestSuite):
143 143 test_support.pickeling(self, gammalib.GModelSpectralGauss())
144 144 test_support.pickeling(self, gammalib.GModelSpectralLogParabola())
145 145 test_support.pickeling(self, gammalib.GModelSpectralMultiplicative())
  146 + test_support.pickeling(self, gammalib.GModelSpectralBins())
146 147 test_support.pickeling(self, gammalib.GModelSpectralNodes())
147 148 test_support.pickeling(self, gammalib.GModelSpectralPlaw())
148 149 test_support.pickeling(self, gammalib.GModelSpectralPlawEnergyFlux())
... ... @@ -169,6 +170,7 @@ class Test(gammalib.GPythonTestSuite):
169 170 cnst = gammalib.GModelTemporalConst(3.0)
170 171 sky = gammalib.GModelSky(ptsrc, plaw, cnst)
171 172 engs = gammalib.GEnergies(10,pivot,gammalib.GEnergy(10.0,'TeV'))
  173 + ebds = gammalib.GEbounds(10,gammalib.GEnergy(1.0,'TeV'),gammalib.GEnergy(10.0,'TeV'))
172 174 circle = gammalib.GSkyRegionCircle(83.6331, 22.0145, 3.0)
173 175 spacom = gammalib.GModelSpatialComposite()
174 176 spacom.append(ptsrc, 'Src1')
... ... @@ -181,6 +183,7 @@ class Test(gammalib.GPythonTestSuite):
181 183 spemul = gammalib.GModelSpectralMultiplicative()
182 184 spemul.append(plaw, 'Src1')
183 185 spemul.append(plaw, 'Src2')
  186 + spebin = gammalib.GModelSpectralBins(plaw, ebds)
184 187 spenod = gammalib.GModelSpectralNodes(plaw, engs)
185 188 spamap = gammalib.GModelSpatialDiffuseMap(cube, 2.0)
186 189 spamap.mc_cone(circle)
... ... @@ -231,6 +234,7 @@ class Test(gammalib.GPythonTestSuite):
231 234 test_support.pickeling(self, gammalib.GModelSpectralGauss(1.0,pivot,pivot))
232 235 test_support.pickeling(self, gammalib.GModelSpectralLogParabola(1.0,-2.0,pivot,0.9))
233 236 test_support.pickeling(self, gammalib.GModelSpectralMultiplicative(spemul))
  237 + test_support.pickeling(self, gammalib.GModelSpectralBins(spebin))
234 238 test_support.pickeling(self, gammalib.GModelSpectralNodes(spenod))
235 239 test_support.pickeling(self, gammalib.GModelSpectralPlaw(plaw))
236 240 test_support.pickeling(self, gammalib.GModelSpectralPlawEnergyFlux(1.0,-2.0,pivot,pivot))
... ...