Olga Sergijenko / gammalib

Commit 52e5c22c1f076fc5bbca078b1f312718265df448

Authored by Alexander Ziegler
1 parent 4bd3d7d4

Implementations for new spectral class GModelSpectralExpPlaw2

Inline Side-by-side
Showing 11 changed files with 1564 additions and 0 deletions
doc/source/users/user_manual/modules/model.rst
  View file @52e5c22
... ... @@ -423,6 +423,33 @@ The XML format for specifying an exponentially cut-off power law is:
423 423 <parameter name="Scale" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
424 424 </spectrum>
425 425  
  426 +An alternative exponentially cut-off power law function is defined by the
  427 +:doxy:`GModelSpectralExpPlaw2` class which makes use of the inverse of the
  428 +cut-off energy for function parametrisation:
  429 +
  430 +.. math::
  431 + \frac{dN}{dE} = k_0 \left( \frac{E}{E_0} \right)^{\gamma}
  432 + \exp \left( -\lambda E \right)
  433 +
  434 +where the parameters in the XML definition have the following mappings:
  435 +
  436 +* :math:`k_0` = ``Prefactor``
  437 +* :math:`\gamma` = ``Index``
  438 +* :math:`E_0` = ``Scale``
  439 +* :math:`\lambda` = ``lambda``
  440 +
  441 +The XML format for specifying an exponentially cut-off power law using this
  442 +alternative parametrisation is:
  443 +
  444 +.. code-block:: xml
  445 +
  446 + <spectrum type="ExpCutoff2">
  447 + <parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
  448 + <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
  449 + <parameter name="lambda" scale="1e-6" value="1.0" min="0.0" max="100.0" free="1"/>
  450 + <parameter name="Scale" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
  451 + </spectrum>
  452 +
426 453  
427 454 Super exponentially cut-off power law
428 455 =====================================
... ...
include/GModelSpectralExpPlaw2.hpp 0 → 100644
  View file @52e5c22
  1 +/***************************************************************************
  2 + * GModelSpectralExpPlaw.hpp - Exponential cut off power law model *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2010-2016 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 GModelSpectralExpPlaw2.hpp
  23 + * @brief Exponential cut off power law spectral class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +
  27 +#ifndef GMODELSPECTRALEXPPLAW2_HPP
  28 +#define GMODELSPECTRALEXPPLAW2_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GModelSpectral.hpp"
  33 +#include "GFunction.hpp"
  34 +#include "GModelPar.hpp"
  35 +#include "GEnergy.hpp"
  36 +
  37 +/* __ Forward declarations _______________________________________________ */
  38 +class GRan;
  39 +class GTime;
  40 +class GXmlElement;
  41 +
  42 +
  43 +/***********************************************************************//**
  44 + * @class GModelSpectralExpPlaw2
  45 + *
  46 + * @brief Exponential cut off power law spectral class
  47 + *
  48 + * This class implements a power law spectrum with exponential cut off. The
  49 + * model is defined by
  50 + *
  51 + * \f[
  52 + * S_{\rm E}(E | t) = {\tt m\_norm}
  53 + * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index}
  54 + * \exp \left( -{\tt m\_lamda}*E \right)
  55 + * \f]
  56 + *
  57 + * where
  58 + * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  59 + * - \f${\tt m\_pivot}\f$ is the pivot energy,
  60 + * - \f${\tt m\_index}\f$ is the spectral index, and
  61 + * - \f${\tt m\_lambda}\f$ is the cut-off parameter (1/cut-off energy).
  62 + ***************************************************************************/
  63 +class GModelSpectralExpPlaw2 : public GModelSpectral {
  64 +
  65 +public:
  66 + // Constructors and destructors
  67 + GModelSpectralExpPlaw2(void);
  68 + explicit GModelSpectralExpPlaw2(const double& prefactor,
  69 + const double& index,
  70 + const GEnergy& pivot,
  71 + const double& lambda);
  72 + explicit GModelSpectralExpPlaw2(const GXmlElement& xml);
  73 + GModelSpectralExpPlaw2(const GModelSpectralExpPlaw2& model);
  74 + virtual ~GModelSpectralExpPlaw2(void);
  75 +
  76 + // Operators
  77 + virtual GModelSpectralExpPlaw2& operator=(const GModelSpectralExpPlaw2& model);
  78 +
  79 + // Implemented pure virtual methods
  80 + virtual void clear(void);
  81 + virtual GModelSpectralExpPlaw2* clone(void) const;
  82 + virtual std::string classname(void) const;
  83 + virtual std::string type(void) const;
  84 + virtual double eval(const GEnergy& srcEng,
  85 + const GTime& srcTime = GTime()) const;
  86 + virtual double eval_gradients(const GEnergy& srcEng,
  87 + const GTime& srcTime = GTime());
  88 + virtual double flux(const GEnergy& emin,
  89 + const GEnergy& emax) const;
  90 + virtual double eflux(const GEnergy& emin,
  91 + const GEnergy& emax) const;
  92 + virtual GEnergy mc(const GEnergy& emin,
  93 + const GEnergy& emax,
  94 + const GTime& time,
  95 + GRan& ran) const;
  96 + virtual void read(const GXmlElement& xml);
  97 + virtual void write(GXmlElement& xml) const;
  98 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  99 +
  100 + // Other methods
  101 + double prefactor(void) const;
  102 + void prefactor(const double& prefactor);
  103 + double index(void) const;
  104 + void index(const double& index);
  105 + double lambda(void) const;
  106 + void lambda(const double& lambda);
  107 + GEnergy pivot(void) const;
  108 + void pivot(const GEnergy& pivot);
  109 +
  110 +protected:
  111 + // Protected methods
  112 + void init_members(void);
  113 + void copy_members(const GModelSpectralExpPlaw2& model);
  114 + void free_members(void);
  115 + void update_eval_cache(const GEnergy& energy) const;
  116 + void update_mc_cache(const GEnergy& emin, const GEnergy& emax) const;
  117 +
  118 + // Photon flux integration kernel
  119 + class flux_kernel : public GFunction {
  120 + public:
  121 + flux_kernel(const double& norm,
  122 + const double& index,
  123 + const double& pivot,
  124 + const double& lambda) :
  125 + m_norm(norm),
  126 + m_index(index),
  127 + m_inv_pivot(1.0/pivot),
  128 + m_lambda(lambda) {}
  129 + double eval(const double& eng);
  130 + protected:
  131 + double m_norm; //!< Normalization
  132 + double m_index; //!< Index
  133 + double m_inv_pivot; //!< 1 / Pivot energy
  134 + double m_lambda; //!< Cut-off parameter
  135 + };
  136 +
  137 + // Energy flux integration kernel
  138 + class eflux_kernel : public GFunction {
  139 + public:
  140 + eflux_kernel(const double& norm,
  141 + const double& index,
  142 + const double& pivot,
  143 + const double& lambda) :
  144 + m_norm(norm),
  145 + m_index(index),
  146 + m_inv_pivot(1.0/pivot),
  147 + m_lambda(lambda) {}
  148 + double eval(const double& eng);
  149 + protected:
  150 + double m_norm; //!< Normalization
  151 + double m_index; //!< Index
  152 + double m_inv_pivot; //!< 1 / Pivot energy
  153 + double m_lambda; //!< Cut-off parameter
  154 + };
  155 +
  156 + // Protected members
  157 + GModelPar m_norm; //!< Normalization factor
  158 + GModelPar m_index; //!< Spectral index
  159 + GModelPar m_lambda; //!< Cut-off parameter
  160 + GModelPar m_pivot; //!< Pivot energy
  161 +
  162 + // Cached members used for pre-computations
  163 + mutable GEnergy m_last_energy; //!< Last energy value
  164 + mutable double m_last_index; //!< Last index parameter
  165 + mutable double m_last_lambda; //!< Last cut-off parameter
  166 + mutable double m_last_pivot; //!< Last pivot parameter
  167 + mutable double m_last_e_norm; //!< Last E/Epivot value
  168 + mutable double m_last_e_lambda; //!< Last E*lambda value
  169 + mutable double m_last_power; //!< Last power value
  170 + mutable double m_mc_emin; //!< Minimum energy
  171 + mutable double m_mc_emax; //!< Maximum energy
  172 + mutable double m_mc_exponent; //!< Exponent (index+1)
  173 + mutable double m_mc_pow_emin; //!< Power of minimum energy
  174 + mutable double m_mc_pow_ewidth; //!< Power of energy width
  175 +};
  176 +
  177 +
  178 +/***********************************************************************//**
  179 + * @brief Return class name
  180 + *
  181 + * @return String containing the class name ("GModelSpectralExpPlaw2").
  182 + ***************************************************************************/
  183 +inline
  184 +std::string GModelSpectralExpPlaw2::classname(void) const
  185 +{
  186 + return ("GModelSpectralExpPlaw2");
  187 +}
  188 +
  189 +
  190 +/***********************************************************************//**
  191 + * @brief Return model type
  192 + *
  193 + * @return "ExpCutoff2".
  194 + *
  195 + * Returns the type of the exponentially cut off power law model.
  196 + ***************************************************************************/
  197 +inline
  198 +std::string GModelSpectralExpPlaw2::type(void) const
  199 +{
  200 + return "ExpCutoff2";
  201 +}
  202 +
  203 +
  204 +/***********************************************************************//**
  205 + * @brief Return pre factor
  206 + *
  207 + * @return Pre factor (ph/cm2/s/MeV).
  208 + *
  209 + * Returns the pre factor.
  210 + ***************************************************************************/
  211 +inline
  212 +double GModelSpectralExpPlaw2::prefactor(void) const
  213 +{
  214 + return (m_norm.value());
  215 +}
  216 +
  217 +
  218 +/***********************************************************************//**
  219 + * @brief Set pre factor
  220 + *
  221 + * @param[in] prefactor Pre factor (ph/cm2/s/MeV).
  222 + *
  223 + * Sets the pre factor.
  224 + ***************************************************************************/
  225 +inline
  226 +void GModelSpectralExpPlaw2::prefactor(const double& prefactor)
  227 +{
  228 + m_norm.value(prefactor);
  229 + return;
  230 +}
  231 +
  232 +
  233 +/***********************************************************************//**
  234 + * @brief Return power law index
  235 + *
  236 + * @return Power law index.
  237 + *
  238 + * Returns the power law index.
  239 + ***************************************************************************/
  240 +inline
  241 +double GModelSpectralExpPlaw2::index(void) const
  242 +{
  243 + return (m_index.value());
  244 +}
  245 +
  246 +
  247 +/***********************************************************************//**
  248 + * @brief Set power law index
  249 + *
  250 + * @param[in] index Power law index.
  251 + *
  252 + * Sets the power law index.
  253 + ***************************************************************************/
  254 +inline
  255 +void GModelSpectralExpPlaw2::index(const double& index)
  256 +{
  257 + m_index.value(index);
  258 + return;
  259 +}
  260 +
  261 +
  262 +/***********************************************************************//**
  263 + * @brief Return pivot energy
  264 + *
  265 + * @return Pivot energy.
  266 + *
  267 + * Returns the pivot energy.
  268 + ***************************************************************************/
  269 +inline
  270 +GEnergy GModelSpectralExpPlaw2::pivot(void) const
  271 +{
  272 + GEnergy energy;
  273 + energy.MeV(m_pivot.value());
  274 + return energy;
  275 +}
  276 +
  277 +
  278 +/***********************************************************************//**
  279 + * @brief Set pivot energy
  280 + *
  281 + * @param[in] pivot Pivot energy.
  282 + *
  283 + * Sets the pivot energy.
  284 + ***************************************************************************/
  285 +inline
  286 +void GModelSpectralExpPlaw2::pivot(const GEnergy& pivot)
  287 +{
  288 + m_pivot.value(pivot.MeV());
  289 + return;
  290 +}
  291 +
  292 +
  293 +/***********************************************************************//**
  294 + * @brief Return exponential cut-off parameter
  295 + *
  296 + * @return Exponential cut-off parameter.
  297 + *
  298 + * Returns the exponential cut-off parameter.
  299 + ***************************************************************************/
  300 +inline
  301 +double GModelSpectralExpPlaw2::lambda(void) const
  302 +{
  303 + return (m_lambda.value());
  304 +}
  305 +
  306 +
  307 +/***********************************************************************//**
  308 + * @brief Set exponential cut-off parameter
  309 + *
  310 + * @param[in] lambda Exponential cut-off parameter.
  311 + *
  312 + * Sets the exponential cut-off parameter.
  313 + ***************************************************************************/
  314 +inline
  315 +void GModelSpectralExpPlaw2::lambda(const double& lambda)
  316 +{
  317 + m_lambda.value(lambda);
  318 + return;
  319 +}
  320 +
  321 +#endif /* GMODELSPECTRALEXPPLAW2_HPP */
... ...
include/GammaLib.hpp
  View file @52e5c22
... ... @@ -209,6 +209,7 @@
209 209 #include "GModelSpectralBrokenPlaw.hpp"
210 210 #include "GModelSpectralConst.hpp"
211 211 #include "GModelSpectralExpPlaw.hpp"
  212 +#include "GModelSpectralExpPlaw2.hpp"
212 213 #include "GModelSpectralSuperExpPlaw.hpp"
213 214 #include "GModelSpectralFunc.hpp"
214 215 #include "GModelSpectralGauss.hpp"
... ...
include/Makefile.am
  View file @52e5c22
... ... @@ -146,6 +146,7 @@ pkginclude_HEADERS = GammaLib.hpp \
146 146 GModelSpectralBrokenPlaw.hpp \
147 147 GModelSpectralConst.hpp \
148 148 GModelSpectralExpPlaw.hpp \
  149 + GModelSpectralExpPlaw2.hpp \
149 150 GModelSpectralSuperExpPlaw.hpp \
150 151 GModelSpectralFunc.hpp \
151 152 GModelSpectralGauss.hpp \
... ...
pyext/GModelSpectralExpPlaw2.i 0 → 100644
  View file @52e5c22
  1 +/***************************************************************************
  2 + * GModelSpectralExpPlaw.i - Exponential cut off power law model *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2009-2016 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 GModelSpectralExpPlaw2.i
  23 + * @brief Exponential cut off power law spectral class interface definition
  24 + * @author Juergen Knoedlseder
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelSpectralExpPlaw2.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GModelSpectralExpPlaw2
  34 + *
  35 + * @brief Exponential cut off power law spectral class
  36 + ***************************************************************************/
  37 +class GModelSpectralExpPlaw2 : public GModelSpectral {
  38 +
  39 +public:
  40 + // Constructors and destructors
  41 + GModelSpectralExpPlaw2(void);
  42 + explicit GModelSpectralExpPlaw2(const double& prefactor,
  43 + const double& index,
  44 + const GEnergy& pivot,
  45 + const double& lambda);
  46 + explicit GModelSpectralExpPlaw2(const GXmlElement& xml);
  47 + GModelSpectralExpPlaw2(const GModelSpectralExpPlaw2& model);
  48 + virtual ~GModelSpectralExpPlaw2(void);
  49 +
  50 + // Operators
  51 + virtual GModelSpectralExpPlaw2& operator=(const GModelSpectralExpPlaw2& model);
  52 +
  53 + // Implemented pure virtual methods
  54 + virtual void clear(void);
  55 + virtual GModelSpectralExpPlaw2* clone(void) const;
  56 + virtual std::string classname(void) const;
  57 + virtual std::string type(void) const;
  58 + virtual double eval(const GEnergy& srcEng,
  59 + const GTime& srcTime = GTime()) const;
  60 + virtual double eval_gradients(const GEnergy& srcEng,
  61 + const GTime& srcTime = GTime());
  62 + virtual double flux(const GEnergy& emin,
  63 + const GEnergy& emax) const;
  64 + virtual double eflux(const GEnergy& emin,
  65 + const GEnergy& emax) const;
  66 + virtual GEnergy mc(const GEnergy& emin,
  67 + const GEnergy& emax,
  68 + const GTime& time,
  69 + GRan& ran) const;
  70 + virtual void read(const GXmlElement& xml);
  71 + virtual void write(GXmlElement& xml) const;
  72 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  73 +
  74 + // Other methods
  75 + double prefactor(void) const;
  76 + void prefactor(const double& prefactor);
  77 + double index(void) const;
  78 + void index(const double& index);
  79 + double lambda(void) const;
  80 + void lambda(const double& lambda);
  81 + GEnergy pivot(void) const;
  82 + void pivot(const GEnergy& pivot);
  83 +};
  84 +
  85 +
  86 +/***********************************************************************//**
  87 + * @brief GModelSpectralExpPlaw2 class extension
  88 + ***************************************************************************/
  89 +%extend GModelSpectralExpPlaw2 {
  90 + GModelSpectralExpPlaw2 copy() {
  91 + return (*self);
  92 + }
  93 +};
... ...
pyext/gammalib/model.i
  View file @52e5c22
... ... @@ -146,6 +146,7 @@
146 146 %include "GModelSpectralBrokenPlaw.i"
147 147 %include "GModelSpectralConst.i"
148 148 %include "GModelSpectralExpPlaw.i"
  149 +%include "GModelSpectralExpPlaw2.i"
149 150 %include "GModelSpectralSuperExpPlaw.i"
150 151 %include "GModelSpectralFunc.i"
151 152 %include "GModelSpectralGauss.i"
... ...
src/model/GModelSpectralExpPlaw2.cpp 0 → 100644
  View file @52e5c22
  1 +/***************************************************************************
  2 + * GModelSpectralExpPlaw.cpp - Exponential cut off power law model *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2010-2016 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 GModelSpectralExpPlaw2.cpp
  23 + * @brief Exponential cut off power law spectral 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 "GIntegral.hpp"
  36 +#include "GModelSpectralExpPlaw2.hpp"
  37 +#include "GModelSpectralRegistry.hpp"
  38 +
  39 +/* __ Constants __________________________________________________________ */
  40 +
  41 +/* __ Globals ____________________________________________________________ */
  42 +const GModelSpectralExpPlaw2 g_spectral_expplaw2_seed;
  43 +const GModelSpectralRegistry g_spectral_expplaw2_registry(&g_spectral_expplaw2_seed);
  44 +
  45 +/* __ Method name definitions ____________________________________________ */
  46 +#define G_MC "GModelSpectralExpPlaw2::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
  47 +#define G_READ "GModelSpectralExpPlaw2::read(GXmlElement&)"
  48 +#define G_WRITE "GModelSpectralExpPlaw2::write(GXmlElement&)"
  49 +
  50 +/* __ Macros _____________________________________________________________ */
  51 +
  52 +/* __ Coding definitions _________________________________________________ */
  53 +
  54 +/* __ Debug definitions __________________________________________________ */
  55 +
  56 +
  57 +/*==========================================================================
  58 + = =
  59 + = Constructors/destructors =
  60 + = =
  61 + ==========================================================================*/
  62 +
  63 +/***********************************************************************//**
  64 + * @brief Void constructor
  65 + ***************************************************************************/
  66 +GModelSpectralExpPlaw2::GModelSpectralExpPlaw2(void) : GModelSpectral()
  67 +{
  68 + // Initialise members
  69 + init_members();
  70 +
  71 + // Return
  72 + return;
  73 +}
  74 +
  75 +
  76 +/***********************************************************************//**
  77 + * @brief Constructor
  78 + *
  79 + * @param[in] prefactor Pre factor normalization (ph/cm2/s/MeV).
  80 + * @param[in] index Power law index.
  81 + * @param[in] pivot Pivot energy.
  82 + * @param[in] lambda Cut-off parameter.
  83 + *
  84 + * Construct an exponentially cut off power law from
  85 + * - a prefactor value (in units of ph/cm2/s/MeV)
  86 + * - a spectral index,
  87 + * - a pivot energy, and
  88 + * - a cut-off parameter lambda.
  89 + ***************************************************************************/
  90 +GModelSpectralExpPlaw2::GModelSpectralExpPlaw2(const double& prefactor,
  91 + const double& index,
  92 + const GEnergy& pivot,
  93 + const double & lambda) :
  94 + GModelSpectral()
  95 +{
  96 + // Initialise members
  97 + init_members();
  98 +
  99 + // Set parameters
  100 + m_norm.value(prefactor);
  101 + m_index.value(index);
  102 + m_pivot.value(pivot.MeV()); // Internally stored in MeV
  103 + m_lambda.value(lambda); // Internally stored in 1/MeV
  104 +
  105 + // Autoscale parameters
  106 + autoscale();
  107 +
  108 + // Return
  109 + return;
  110 +}
  111 +
  112 +
  113 +/***********************************************************************//**
  114 + * @brief XML constructor
  115 + *
  116 + * @param[in] xml XML element.
  117 + *
  118 + * Constructs an exponentially cut off power law spectral model by extracting
  119 + * information from an XML element. See the read() method for more
  120 + * information about the expected structure of the XML element.
  121 + ***************************************************************************/
  122 +GModelSpectralExpPlaw2::GModelSpectralExpPlaw2(const GXmlElement& xml) :
  123 + GModelSpectral()
  124 +{
  125 + // Initialise members
  126 + init_members();
  127 +
  128 + // Read information from XML element
  129 + read(xml);
  130 +
  131 + // Return
  132 + return;
  133 +}
  134 +
  135 +
  136 +/***********************************************************************//**
  137 + * @brief Copy constructor
  138 + *
  139 + * @param[in] model Exponentially cut off power law model.
  140 + ***************************************************************************/
  141 +GModelSpectralExpPlaw2::GModelSpectralExpPlaw2(const GModelSpectralExpPlaw2& model) :
  142 + GModelSpectral(model)
  143 +{
  144 + // Initialise members
  145 + init_members();
  146 +
  147 + // Copy members
  148 + copy_members(model);
  149 +
  150 + // Return
  151 + return;
  152 +}
  153 +
  154 +
  155 +/***********************************************************************//**
  156 + * @brief Destructor
  157 + ***************************************************************************/
  158 +GModelSpectralExpPlaw2::~GModelSpectralExpPlaw2(void)
  159 +{
  160 + // Free members
  161 + free_members();
  162 +
  163 + // Return
  164 + return;
  165 +}
  166 +
  167 +
  168 +/*==========================================================================
  169 + = =
  170 + = Operators =
  171 + = =
  172 + ==========================================================================*/
  173 +
  174 +/***********************************************************************//**
  175 + * @brief Assignment operator
  176 + *
  177 + * @param[in] model Exponentially cut off power law model.
  178 + * @return Exponentially cut off power law model.
  179 + ***************************************************************************/
  180 +GModelSpectralExpPlaw2& GModelSpectralExpPlaw2::operator=(const GModelSpectralExpPlaw2& model)
  181 +{
  182 + // Execute only if object is not identical
  183 + if (this != &model) {
  184 +
  185 + // Copy base class members
  186 + this->GModelSpectral::operator=(model);
  187 +
  188 + // Free members
  189 + free_members();
  190 +
  191 + // Initialise members
  192 + init_members();
  193 +
  194 + // Copy members
  195 + copy_members(model);
  196 +
  197 + } // endif: object was not identical
  198 +
  199 + // Return
  200 + return *this;
  201 +}
  202 +
  203 +
  204 +/*==========================================================================
  205 + = =
  206 + = Public methods =
  207 + = =
  208 + ==========================================================================*/
  209 +
  210 +/***********************************************************************//**
  211 + * @brief Clear exponentially cut off power law model
  212 + ***************************************************************************/
  213 +void GModelSpectralExpPlaw2::clear(void)
  214 +{
  215 + // Free class members (base and derived classes, derived class first)
  216 + free_members();
  217 + this->GModelSpectral::free_members();
  218 +
  219 + // Initialise members
  220 + this->GModelSpectral::init_members();
  221 + init_members();
  222 +
  223 + // Return
  224 + return;
  225 +}
  226 +
  227 +
  228 +/***********************************************************************//**
  229 + * @brief Clone exponentially cut off power law model
  230 + *
  231 + * @return Pointer to deep copy of exponentially cut off power law model.
  232 + ***************************************************************************/
  233 +GModelSpectralExpPlaw2* GModelSpectralExpPlaw2::clone(void) const
  234 +{
  235 + // Clone exponentially cut off power law model
  236 + return new GModelSpectralExpPlaw2(*this);
  237 +}
  238 +
  239 +
  240 +/***********************************************************************//**
  241 + * @brief Evaluate function
  242 + *
  243 + * @param[in] srcEng True photon energy.
  244 + * @param[in] srcTime True photon arrival time.
  245 + * @return Model value.
  246 + *
  247 + * Evaluates
  248 + *
  249 + * \f[
  250 + * S_{\rm E}(E | t) = {\tt m\_norm}
  251 + * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index}
  252 + * \exp \left( -{\tt m\_lamda}*E \right)
  253 + * \f]
  254 + *
  255 + * where
  256 + * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  257 + * - \f${\tt m\_pivot}\f$ is the pivot energy,
  258 + * - \f${\tt m\_index}\f$ is the spectral index, and
  259 + * - \f${\tt m\_lambda}\f$ is the cut-off parameter.
  260 + *
  261 + * @todo The method expects that pivot!=0. Otherwise Inf or NaN
  262 + * may result. We should add a test that prevents using invalid
  263 + * values.
  264 + ***************************************************************************/
  265 +double GModelSpectralExpPlaw2::eval(const GEnergy& srcEng,
  266 + const GTime& srcTime) const
  267 +{
  268 + // Update the evaluation cache
  269 + update_eval_cache(srcEng);
  270 +
  271 + // Compute function value
  272 + double value = m_norm.value() * m_last_power;
  273 +
  274 + // Compile option: Check for NaN/Inf
  275 + #if defined(G_NAN_CHECK)
  276 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  277 + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval";
  278 + std::cout << "(srcEng=" << srcEng;
  279 + std::cout << ", srcTime=" << srcTime << "):";
  280 + std::cout << " NaN/Inf encountered";
  281 + std::cout << " (value=" << value;
  282 + std::cout << ", power=" << m_last_power;
  283 + std::cout << ")" << std::endl;
  284 + }
  285 + #endif
  286 +
  287 + // Return
  288 + return value;
  289 +}
  290 +
  291 +
  292 +/***********************************************************************//**
  293 + * @brief Evaluate function and gradients
  294 + *
  295 + * @param[in] srcEng True photon energy.
  296 + * @param[in] srcTime True photon arrival time.
  297 + * @return Model value.
  298 + *
  299 + * Evaluates
  300 + *
  301 + * \f[
  302 + * S_{\rm E}(E | t) = {\tt m\_norm}
  303 + * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index}
  304 + * \exp \left( -{\tt m\_lamda}*E \right)
  305 + * \f]
  306 + *
  307 + * where
  308 + * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  309 + * - \f${\tt m\_pivot}\f$ is the pivot energy,
  310 + * - \f${\tt m\_index}\f$ is the spectral index, and
  311 + * - \f${\tt m\_lambda}\f$ is the cut-off parameter.
  312 + *
  313 + * The method also evaluates the partial derivatives of the model with
  314 + * respect to the parameters using
  315 + *
  316 + * \f[
  317 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} =
  318 + * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}}
  319 + * \f]
  320 + *
  321 + * \f[
  322 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} =
  323 + * S_{\rm E}(E | t) \, \ln(E/{\tt m_pivot})
  324 + * \f]
  325 + *
  326 + * \f[
  327 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_lambda}} =
  328 + * S_{\rm E}(E | t) \, \left( -E \right)
  329 + * \f]
  330 + *
  331 + * \f[
  332 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_pivot}} =
  333 + * -S_{\rm E}(E | t) \,
  334 + * \left( \frac{{\tt m\_index}}{{\tt m\_pivot}} \right)
  335 + * \f]
  336 + *
  337 + * @todo The method expects that pivot!=0. Otherwise Inf or NaN
  338 + * may result. We should add a test that prevents using invalid
  339 + * values.
  340 + ***************************************************************************/
  341 +double GModelSpectralExpPlaw2::eval_gradients(const GEnergy& srcEng,
  342 + const GTime& srcTime)
  343 +{
  344 + // Update the evaluation cache
  345 + update_eval_cache(srcEng);
  346 +
  347 + // Compute function value
  348 + double value = m_norm.value() * m_last_power;
  349 +
  350 + // Compute partial derivatives with respect to the parameter factor
  351 + // values. The partial derivatives with respect to the parameter
  352 + // values are obtained by division by the scale factor.
  353 + double g_norm = (m_norm.is_free())
  354 + ? m_norm.scale() * m_last_power : 0.0;
  355 + double g_index = (m_index.is_free())
  356 + ? value * m_index.scale() * std::log(m_last_e_norm) : 0.0;
  357 + double g_lambda = (m_lambda.is_free())
  358 + ? -value * m_lambda.scale() * m_last_energy.MeV() : 0.0;
  359 + double g_pivot = (m_pivot.is_free())
  360 + ? -value * m_last_index / m_pivot.factor_value() : 0.0;
  361 +
  362 + // Set gradients
  363 + m_norm.factor_gradient(g_norm);
  364 + m_index.factor_gradient(g_index);
  365 + m_lambda.factor_gradient(g_lambda);
  366 + m_pivot.factor_gradient(g_pivot);
  367 +
  368 + // Compile option: Check for NaN/Inf
  369 + #if defined(G_NAN_CHECK)
  370 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
  371 + std::cout << "*** ERROR: GModelSpectralExpPlaw::eval_gradients";
  372 + std::cout << "(srcEng=" << srcEng;
  373 + std::cout << ", srcTime=" << srcTime << "):";
  374 + std::cout << " NaN/Inf encountered";
  375 + std::cout << " (value=" << value;
  376 + std::cout << ", e_norm=" << m_last_e_norm;
  377 + std::cout << ", e_lambda=" << m_last_e_lambda;
  378 + std::cout << ", power=" << m_last_power;
  379 + std::cout << ")" << std::endl;
  380 + }
  381 + #endif
  382 +
  383 + // Return
  384 + return value;
  385 +}
  386 +
  387 +
  388 +/***********************************************************************//**
  389 + * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
  390 + *
  391 + * @param[in] emin Minimum photon energy.
  392 + * @param[in] emax Maximum photon energy.
  393 + * @return Photon flux (ph/cm2/s).
  394 + *
  395 + * Computes
  396 + *
  397 + * \f[
  398 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
  399 + * \f]
  400 + *
  401 + * where
  402 + * - [@p emin, @p emax] is an energy interval, and
  403 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  404 + * The integration is done numerically.
  405 + ***************************************************************************/
  406 +double GModelSpectralExpPlaw2::flux(const GEnergy& emin,
  407 + const GEnergy& emax) const
  408 +{
  409 + // Initialise flux
  410 + double flux = 0.0;
  411 +
  412 + // Compute only if integration range is valid
  413 + if (emin < emax) {
  414 +
  415 + // Setup integration kernel
  416 + flux_kernel integrand(m_norm.value(), m_index.value(),
  417 + m_pivot.value(), m_lambda.value());
  418 + GIntegral integral(&integrand);
  419 +
  420 + // Get integration boundaries in MeV
  421 + double e_min = emin.MeV();
  422 + double e_max = emax.MeV();
  423 +
  424 + // Perform integration
  425 + flux = integral.romberg(e_min, e_max);
  426 +
  427 + } // endif: integration range was valid
  428 +
  429 + // Return
  430 + return flux;
  431 +}
  432 +
  433 +
  434 +/***********************************************************************//**
  435 + * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
  436 + *
  437 + * @param[in] emin Minimum photon energy.
  438 + * @param[in] emax Maximum photon energy.
  439 + * @return Energy flux (erg/cm2/s).
  440 + *
  441 + * Computes
  442 + *
  443 + * \f[
  444 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
  445 + * \f]
  446 + *
  447 + * where
  448 + * - [@p emin, @p emax] is an energy interval, and
  449 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  450 + * The integration is done numerically.
  451 + ***************************************************************************/
  452 +double GModelSpectralExpPlaw2::eflux(const GEnergy& emin,
  453 + const GEnergy& emax) const
  454 +{
  455 + // Initialise flux
  456 + double eflux = 0.0;
  457 +
  458 + // Compute only if integration range is valid
  459 + if (emin < emax) {
  460 +
  461 + // Setup integration kernel
  462 + eflux_kernel integrand(m_norm.value(), m_index.value(),
  463 + m_pivot.value(), m_lambda.value());
  464 + GIntegral integral(&integrand);
  465 +
  466 + // Get integration boundaries in MeV
  467 + double e_min = emin.MeV();
  468 + double e_max = emax.MeV();
  469 +
  470 + // Perform integration
  471 + eflux = integral.romberg(e_min, e_max);
  472 +
  473 + // Convert from MeV/cm2/s to erg/cm2/s
  474 + eflux *= gammalib::MeV2erg;
  475 +
  476 + } // endif: integration range was valid
  477 +
  478 + // Return
  479 + return eflux;
  480 +}
  481 +
  482 +
  483 +/***********************************************************************//**
  484 + * @brief Returns MC energy between [emin, emax]
  485 + *
  486 + * @param[in] emin Minimum photon energy.
  487 + * @param[in] emax Maximum photon energy.
  488 + * @param[in] time True photon arrival time.
  489 + * @param[in,out] ran Random number generator.
  490 + * @return Energy.
  491 + *
  492 + * @exception GException::erange_invalid
  493 + * Energy range is invalid (emin < emax required).
  494 + *
  495 + * Simulates a random energy in the interval [emin, emax] for an
  496 + * exponentially cut off power law. The simulation is done using a rejection
  497 + * method. First, a random energy within [emin, emax] is drawn from an
  498 + * power law distribution. Then the energy is accepted or rejected based
  499 + * on an acceptance fraction that is computed from the exponential cut off.
  500 + ***************************************************************************/
  501 +GEnergy GModelSpectralExpPlaw2::mc(const GEnergy& emin,
  502 + const GEnergy& emax,
  503 + const GTime& time,
  504 + GRan& ran) const
  505 +{
  506 + // Throw an exception if energy range is invalid
  507 + if (emin >= emax) {
  508 + throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
  509 + "Minimum energy < maximum energy required.");
  510 + }
  511 +
  512 + // Allocate energy
  513 + GEnergy energy;
  514 +
  515 + // Update cache
  516 + update_mc_cache(emin, emax);
  517 +
  518 + // Initialise energy
  519 + double eng;
  520 +
  521 + // Initialse acceptance fraction
  522 + double acceptance_fraction;
  523 + double inv_ecut = m_lambda.value();
  524 +
  525 + // Use rejection method to draw a random energy. We first draw
  526 + // analytically from a power law, and then compare the power law
  527 + // at the drawn energy to the exponentially cut off function. This
  528 + // gives an acceptance fraction, and we accept the energy only if
  529 + // a uniform random number is <= the acceptance fraction.
  530 + do {
  531 +
  532 + // Get uniform random number
  533 + double u = ran.uniform();
  534 +
  535 + // Case A: Index is not -1
  536 + if (index() != -1.0) {
  537 + if (u > 0.0) {
  538 + eng = std::exp(std::log(u * m_mc_pow_ewidth + m_mc_pow_emin) /
  539 + m_mc_exponent);
  540 + }
  541 + else {
  542 + eng = 0.0;
  543 + }
  544 + }
  545 +
  546 + // Case B: Index is -1
  547 + else {
  548 + eng = std::exp(u * m_mc_pow_ewidth + m_mc_pow_emin);
  549 + }
  550 +
  551 + // Compute acceptance fraction
  552 + acceptance_fraction = std::exp(-eng * inv_ecut);
  553 +
  554 + } while (ran.uniform() > acceptance_fraction);
  555 +
  556 + // Set energy
  557 + energy.MeV(eng);
  558 +
  559 + // Return energy
  560 + return energy;
  561 +}
  562 +
  563 +
  564 +/***********************************************************************//**
  565 + * @brief Read model from XML element
  566 + *
  567 + * @param[in] xml XML element containing power law model information.
  568 + *
  569 + * @exception GException::model_invalid_parnum
  570 + * Invalid number of model parameters found in XML element.
  571 + * @exception GException::model_invalid_parnames
  572 + * Invalid model parameter names found in XML element.
  573 + *
  574 + * Reads the spectral information from an XML element. The format of the XML
  575 + * elements is
  576 + *
  577 + * <spectrum type="ExpCutoff2">
  578 + * <parameter name="Prefactor" scale=".." value=".." min=".." max=".." free=".."/>
  579 + * <parameter name="Index" scale=".." value=".." min=".." max=".." free=".."/>
  580 + * <parameter name="lambda" scale=".." value=".." min=".." max=".." free=".."/>
  581 + * <parameter name="Scale" scale=".." value=".." min=".." max=".." free=".."/>
  582 + * </spectrum>
  583 + ***************************************************************************/
  584 +void GModelSpectralExpPlaw2::read(const GXmlElement& xml)
  585 +{
  586 + // Verify that XML element has exactly 4 parameters
  587 + if (xml.elements() != 4 || xml.elements("parameter") != 4) {
  588 + throw GException::model_invalid_parnum(G_READ, xml,
  589 + "Power law model requires exactly 4 parameters.");
  590 + }
  591 +
  592 + // Extract model parameters
  593 + int npar[] = {0, 0, 0, 0};
  594 + for (int i = 0; i < 4; ++i) {
  595 +
  596 + // Get parameter element
  597 + const GXmlElement* par = xml.element("parameter", i);
  598 +
  599 + // Handle prefactor
  600 + if (par->attribute("name") == "Prefactor") {
  601 + m_norm.read(*par);
  602 + npar[0]++;
  603 + }
  604 +
  605 + // Handle index
  606 + else if (par->attribute("name") == "Index") {
  607 + m_index.read(*par);
  608 + npar[1]++;
  609 + }
  610 +
  611 + // Handle cutoff parameter
  612 + else if (par->attribute("name") == "lambda") {
  613 + m_lambda.read(*par);
  614 + npar[2]++;
  615 + }
  616 +
  617 + // Handle pivot energy
  618 + else if (par->attribute("name") == "Scale") {
  619 + m_pivot.read(*par);
  620 + npar[3]++;
  621 + }
  622 +
  623 + } // endfor: looped over all parameters
  624 +
  625 + // Verify that all parameters were found
  626 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) {
  627 + throw GException::model_invalid_parnames(G_READ, xml,
  628 + "Require \"Prefactor\", \"Index\", \"lambda\" and \"Scale\""
  629 + " parameters.");
  630 + }
  631 +
  632 + // Return
  633 + return;
  634 +}
  635 +
  636 +
  637 +/***********************************************************************//**
  638 + * @brief Write model into XML element
  639 + *
  640 + * @param[in] xml XML element into which model information is written.
  641 + *
  642 + * @exception GException::model_invalid_spectral
  643 + * Existing XML element is not of type "ExpCutoff"
  644 + * @exception GException::model_invalid_parnum
  645 + * Invalid number of model parameters or nodes found in XML element.
  646 + * @exception GException::model_invalid_parnames
  647 + * Invalid model parameter names found in XML element.
  648 + *
  649 + * Writes the spectral information into an XML element. The format of the XML
  650 + * element is
  651 + *
  652 + * <spectrum type="ExpCutoff2">
  653 + * <parameter name="Prefactor" scale=".." value=".." min=".." max=".." free=".."/>
  654 + * <parameter name="Index" scale=".." value=".." min=".." max=".." free=".."/>
  655 + * <parameter name="lambda" scale=".." value=".." min=".." max=".." free=".."/>
  656 + * <parameter name="Scale" scale=".." value=".." min=".." max=".." free=".."/>
  657 + * </spectrum>
  658 + ***************************************************************************/
  659 +void GModelSpectralExpPlaw2::write(GXmlElement& xml) const
  660 +{
  661 + // Set model type
  662 + if (xml.attribute("type") == "") {
  663 + xml.attribute("type", "ExpCutoff2");
  664 + }
  665 +
  666 + // Verify model type
  667 + if (xml.attribute("type") != "ExpCutoff2") {
  668 + throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"),
  669 + "Spectral model is not of type \"ExpCutoff2\".");
  670 + }
  671 +
  672 + // If XML element has 0 nodes then append 4 parameter nodes
  673 + if (xml.elements() == 0) {
  674 + xml.append(GXmlElement("parameter name=\"Prefactor\""));
  675 + xml.append(GXmlElement("parameter name=\"Index\""));
  676 + xml.append(GXmlElement("parameter name=\"lambda\""));
  677 + xml.append(GXmlElement("parameter name=\"Scale\""));
  678 + }
  679 +
  680 + // Verify that XML element has exactly 4 parameters
  681 + if (xml.elements() != 4 || xml.elements("parameter") != 4) {
  682 + throw GException::model_invalid_parnum(G_WRITE, xml,
  683 + "Power law model requires exactly 4 parameters.");
  684 + }
  685 +
  686 + // Set or update model parameter attributes
  687 + int npar[] = {0, 0, 0, 0};
  688 + for (int i = 0; i < 4; ++i) {
  689 +
  690 + // Get parameter element
  691 + GXmlElement* par = xml.element("parameter", i);
  692 +
  693 + // Handle prefactor
  694 + if (par->attribute("name") == "Prefactor") {
  695 + npar[0]++;
  696 + m_norm.write(*par);
  697 + }
  698 +
  699 + // Handle index
  700 + else if (par->attribute("name") == "Index") {
  701 + npar[1]++;
  702 + m_index.write(*par);
  703 + }
  704 +
  705 + // Handle cutoff parameter
  706 + else if (par->attribute("name") == "lambda") {
  707 + npar[2]++;
  708 + m_lambda.write(*par);
  709 + }
  710 +
  711 + // Handle pivot energy
  712 + else if (par->attribute("name") == "Scale") {
  713 + m_pivot.write(*par);
  714 + npar[3]++;
  715 + }
  716 +
  717 + } // endfor: looped over all parameters
  718 +
  719 + // Check of all required parameters are present
  720 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) {
  721 + throw GException::model_invalid_parnames(G_WRITE, xml,
  722 + "Require \"Prefactor\", \"Index\", \"lambda\" and \"Scale\""
  723 + " parameters.");
  724 + }
  725 +
  726 + // Return
  727 + return;
  728 +}
  729 +
  730 +
  731 +/***********************************************************************//**
  732 + * @brief Print model information
  733 + *
  734 + * @param[in] chatter Chattiness (defaults to NORMAL).
  735 + * @return String containing model information.
  736 + ***************************************************************************/
  737 +std::string GModelSpectralExpPlaw2::print(const GChatter& chatter) const
  738 +{
  739 + // Initialise result string
  740 + std::string result;
  741 +
  742 + // Continue only if chatter is not silent
  743 + if (chatter != SILENT) {
  744 +
  745 + // Append header
  746 + result.append("=== GModelSpectralExpPlaw2 ===");
  747 +
  748 + // Append information
  749 + result.append("\n"+gammalib::parformat("Number of parameters"));
  750 + result.append(gammalib::str(size()));
  751 + for (int i = 0; i < size(); ++i) {
  752 + result.append("\n"+m_pars[i]->print(chatter));
  753 + }
  754 +
  755 + } // endif: chatter was not silent
  756 +
  757 + // Return result
  758 + return result;
  759 +}
  760 +
  761 +
  762 +/*==========================================================================
  763 + = =
  764 + = Private methods =
  765 + = =
  766 + ==========================================================================*/
  767 +
  768 +/***********************************************************************//**
  769 + * @brief Initialise class members
  770 + ***************************************************************************/
  771 +void GModelSpectralExpPlaw2::init_members(void)
  772 +{
  773 + // Initialise powerlaw normalisation
  774 + m_norm.clear();
  775 + m_norm.name("Prefactor");
  776 + m_norm.unit("ph/cm2/s/MeV");
  777 + m_norm.scale(1.0);
  778 + m_norm.value(1.0); // default: 1.0
  779 + m_norm.min(0.0); // min: 0.0
  780 + m_norm.free();
  781 + m_norm.gradient(0.0);
  782 + m_norm.has_grad(true);
  783 +
  784 + // Initialise powerlaw index
  785 + m_index.clear();
  786 + m_index.name("Index");
  787 + m_index.scale(1.0);
  788 + m_index.value(-2.0); // default: -2.0
  789 + m_index.range(-10.0,+10.0); // range: [-10,+10]
  790 + m_index.free();
  791 + m_index.gradient(0.0);
  792 + m_index.has_grad(true);
  793 +
  794 + // Initialise cut off parameter
  795 + m_lambda.clear();
  796 + m_lambda.name("lambda");
  797 + m_lambda.unit("1/MeV");
  798 + m_lambda.scale(1.0);
  799 + m_lambda.value(0.001); // default: 0.001
  800 + m_lambda.min(0.0); // min: 0.0
  801 + m_lambda.max(10.0); // max: 10.0
  802 + m_lambda.free();
  803 + m_lambda.gradient(0.0);
  804 + m_lambda.has_grad(true);
  805 +
  806 + // Initialise pivot energy
  807 + m_pivot.clear();
  808 + m_pivot.name("PivotEnergy");
  809 + m_pivot.unit("MeV");
  810 + m_pivot.scale(1.0);
  811 + m_pivot.value(100.0);
  812 + m_pivot.fix();
  813 + m_pivot.gradient(0.0);
  814 + m_pivot.has_grad(true);
  815 +
  816 + // Set parameter pointer(s)
  817 + m_pars.clear();
  818 + m_pars.push_back(&m_norm);
  819 + m_pars.push_back(&m_index);
  820 + m_pars.push_back(&m_lambda);
  821 + m_pars.push_back(&m_pivot);
  822 +
  823 + // Initialise eval cache
  824 + m_last_energy.clear();
  825 + m_last_index = 1.0e30;
  826 + m_last_lambda = 1.0e30;
  827 + m_last_pivot = 1.0e30;
  828 + m_last_e_norm = 0.0;
  829 + m_last_e_lambda = 0.0;
  830 + m_last_power = 0.0;
  831 +
  832 + // Initialise MC cache
  833 + m_mc_emin = 0.0;
  834 + m_mc_emax = 0.0;
  835 + m_mc_exponent = 0.0;
  836 + m_mc_pow_emin = 0.0;
  837 + m_mc_pow_ewidth = 0.0;
  838 +
  839 + // Return
  840 + return;
  841 +}
  842 +
  843 +
  844 +/***********************************************************************//**
  845 + * @brief Copy class members
  846 + *
  847 + * @param[in] model Exponential cut off power law model.
  848 + ***************************************************************************/
  849 +void GModelSpectralExpPlaw2::copy_members(const GModelSpectralExpPlaw2& model)
  850 +{
  851 + // Copy members
  852 + m_norm = model.m_norm;
  853 + m_index = model.m_index;
  854 + m_lambda = model.m_lambda;
  855 + m_pivot = model.m_pivot;
  856 +
  857 + // Set parameter pointer(s)
  858 + m_pars.clear();
  859 + m_pars.push_back(&m_norm);
  860 + m_pars.push_back(&m_index);
  861 + m_pars.push_back(&m_lambda);
  862 + m_pars.push_back(&m_pivot);
  863 +
  864 + // Copy eval cache
  865 + m_last_energy = model.m_last_energy;
  866 + m_last_index = model.m_last_index;
  867 + m_last_lambda = model.m_last_lambda;
  868 + m_last_pivot = model.m_last_pivot;
  869 + m_last_e_norm = model.m_last_e_norm;
  870 + m_last_e_lambda = model.m_last_e_lambda;
  871 + m_last_power = model.m_last_power;
  872 +
  873 + // Copy MC cache
  874 + m_mc_emin = model.m_mc_emin;
  875 + m_mc_emax = model.m_mc_emax;
  876 + m_mc_exponent = model.m_mc_exponent;
  877 + m_mc_pow_emin = model.m_mc_pow_emin;
  878 + m_mc_pow_ewidth = model.m_mc_pow_ewidth;
  879 +
  880 + // Return
  881 + return;
  882 +}
  883 +
  884 +
  885 +/***********************************************************************//**
  886 + * @brief Delete class members
  887 + ***************************************************************************/
  888 +void GModelSpectralExpPlaw2::free_members(void)
  889 +{
  890 + // Return
  891 + return;
  892 +}
  893 +
  894 +
  895 +/***********************************************************************//**
  896 + * @brief Update eval precomputation cache
  897 + *
  898 + * @param[in] energy Energy.
  899 + *
  900 + * Updates the precomputation cache for eval() and eval_gradients() methods.
  901 + ***************************************************************************/
  902 +void GModelSpectralExpPlaw2::update_eval_cache(const GEnergy& energy) const
  903 +{
  904 + // Get parameter values (takes 3 multiplications which are difficult
  905 + // to avoid)
  906 + double index = m_index.value();
  907 + double lambda = m_lambda.value();
  908 + double pivot = m_pivot.value();
  909 +
  910 + // If the energy or one of the parameters index, lambda or pivot
  911 + // energy has changed then recompute the cache
  912 + if ((m_last_energy != energy) ||
  913 + (m_last_index != index) ||
  914 + (m_last_lambda != lambda) ||
  915 + (m_last_pivot != pivot)) {
  916 +
  917 + // Store actual energy and parameter values
  918 + m_last_energy = energy;
  919 + m_last_index = index;
  920 + m_last_lambda = lambda;
  921 + m_last_pivot = pivot;
  922 +
  923 + // Compute and store value
  924 + double eng = energy.MeV();
  925 + m_last_e_norm = eng / m_last_pivot;
  926 + m_last_e_lambda = eng * m_last_lambda;
  927 + m_last_power = std::pow(m_last_e_norm, m_last_index) *
  928 + std::exp(-m_last_e_lambda);
  929 + } // endif: recomputation was required
  930 +
  931 + // Return
  932 + return;
  933 +}
  934 +
  935 +
  936 +/***********************************************************************//**
  937 + * @brief Update Monte Carlo pre computation cache
  938 + *
  939 + * @param[in] emin Minimum photon energy.
  940 + * @param[in] emax Maximum photon energy.
  941 + *
  942 + * Updates the precomputation cache for Monte Carlo simulations.
  943 + ***************************************************************************/
  944 +void GModelSpectralExpPlaw2::update_mc_cache(const GEnergy& emin,
  945 + const GEnergy& emax) const
  946 +
  947 +{
  948 + // Case A: Index is not -1
  949 + if (index() != -1.0) {
  950 +
  951 + // Change in energy boundaries?
  952 + if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) {
  953 + m_mc_emin = emin.MeV();
  954 + m_mc_emax = emax.MeV();
  955 + m_mc_exponent = index() + 1.0;
  956 + m_mc_pow_emin = std::pow(m_mc_emin, m_mc_exponent);
  957 + m_mc_pow_ewidth = std::pow(m_mc_emax, m_mc_exponent) - m_mc_pow_emin;
  958 + }
  959 +
  960 + }
  961 +
  962 + // Case B: Index is -1
  963 + else {
  964 +
  965 + // Change in energy boundaries?
  966 + if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) {
  967 + m_mc_emin = emin.MeV();
  968 + m_mc_emax = emax.MeV();
  969 + m_mc_exponent = 0.0;
  970 + m_mc_pow_emin = std::log(m_mc_emin);
  971 + m_mc_pow_ewidth = std::log(m_mc_emax) - m_mc_pow_emin;
  972 + }
  973 +
  974 + }
  975 +
  976 + // Return
  977 + return;
  978 +}
  979 +
  980 +
  981 +/***********************************************************************//**
  982 + * @brief Kernel for photon flux integration
  983 + *
  984 + * @param[in] energy Energy (MeV).
  985 + ***************************************************************************/
  986 +double GModelSpectralExpPlaw2::flux_kernel::eval(const double& energy)
  987 +{
  988 + // Evaluate function value
  989 + double e_norm = energy * m_inv_pivot;
  990 + double e_cut = energy * m_lambda;
  991 + double power = std::pow(e_norm, m_index) * std::exp(-e_cut);
  992 + double value = m_norm * power;
  993 +
  994 + // Return value
  995 + return value;
  996 +}
  997 +
  998 +
  999 +/***********************************************************************//**
  1000 + * @brief Kernel for energy flux integration
  1001 + *
  1002 + * @param[in] energy Energy (MeV).
  1003 + ***************************************************************************/
  1004 +double GModelSpectralExpPlaw2::eflux_kernel::eval(const double& energy)
  1005 +{
  1006 + // Evaluate function value
  1007 + double e_norm = energy * m_inv_pivot;
  1008 + double e_cut = energy * m_lambda;
  1009 + double power = std::pow(e_norm, m_index) * std::exp(-e_cut);
  1010 + double value = m_norm * power * energy;
  1011 +
  1012 + // Return value
  1013 + return value;
  1014 +}
... ...
src/model/Makefile.am
  View file @52e5c22
... ... @@ -31,6 +31,7 @@ sources = GModelPar.cpp \
31 31 GModelSpectralBrokenPlaw.cpp \
32 32 GModelSpectralConst.cpp \
33 33 GModelSpectralExpPlaw.cpp \
  34 + GModelSpectralExpPlaw2.cpp \
34 35 GModelSpectralSuperExpPlaw.cpp \
35 36 GModelSpectralFunc.cpp \
36 37 GModelSpectralGauss.cpp \
... ...
test/data/model_point_eplaw2.xml 0 → 100644
  View file @52e5c22
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Crab" type="PointSource">
  4 + <spectrum type="ExpCutoff">
  5 + <parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
  6 + <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
  7 + <parameter name="lambda" scale="1e-6" value="1.0" min="0.0" max="100.0" free="1"/>
  8 + <parameter name="Scale" scale="1e6" value="0.3" min="0.01" max="1000.0" free="0"/>
  9 + </spectrum>
  10 + <spatialModel type="SkyDirFunction">
  11 + <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
  12 + <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
  13 + </spatialModel>
  14 + </source>
  15 +</source_library>
... ...
test/test_GModel.cpp
  View file @52e5c22
... ... @@ -53,6 +53,7 @@ void TestGModel::set(void)
53 53 m_xml_model_point_plaw = "data/model_point_plaw.xml";
54 54 m_xml_model_point_plaw2 = "data/model_point_plaw2.xml";
55 55 m_xml_model_point_eplaw = "data/model_point_eplaw.xml";
  56 + m_xml_model_point_eplaw2 = "data/model_point_eplaw2.xml";
56 57 m_xml_model_point_bplaw = "data/model_point_bplaw.xml";
57 58 m_xml_model_point_supeplaw = "data/model_point_supeplaw.xml";
58 59 m_xml_model_point_logparabola = "data/model_point_logparabola.xml";
... ... @@ -89,6 +90,7 @@ void TestGModel::set(void)
89 90 append(static_cast<pfunction>(&TestGModel::test_plaw), "Test GModelSpectralPlaw");
90 91 append(static_cast<pfunction>(&TestGModel::test_plaw2), "Test GModelSpectralPlaw2");
91 92 append(static_cast<pfunction>(&TestGModel::test_eplaw), "Test GModelSpectralExpPlaw");
  93 + append(static_cast<pfunction>(&TestGModel::test_eplaw2), "Test GModelSpectralExpPlaw2");
92 94 append(static_cast<pfunction>(&TestGModel::test_supeplaw), "Test GModelSpectralSuperExpPlaw");
93 95 append(static_cast<pfunction>(&TestGModel::test_bplaw), "Test GModelSpectralBrokenPlaw");
94 96 append(static_cast<pfunction>(&TestGModel::test_logparabola), "Test GModelSpectralLogParabola");
... ... @@ -1731,6 +1733,92 @@ void TestGModel::test_eplaw(void)
1731 1733  
1732 1734  
1733 1735 /***********************************************************************//**
  1736 + * @brief Test GModelSpectralExpPlaw2 class
  1737 + ***************************************************************************/
  1738 +void TestGModel::test_eplaw2(void)
  1739 +{
  1740 + // Test void constructor
  1741 + test_try("Test void constructor");
  1742 + try {
  1743 + GModelSpectralExpPlaw2 model;
  1744 + test_assert(model.type() == "ExpCutoff2",
  1745 + "Model type \"ExpCutoff2\" expected.");
  1746 + test_try_success();
  1747 + }
  1748 + catch (std::exception &e) {
  1749 + test_try_failure(e);
  1750 + }
  1751 +
  1752 + // Test value constructor
  1753 + test_try("Test value constructor");
  1754 + try {
  1755 + GModelSpectralExpPlaw2 model(2.0, -2.1, GEnergy(100.0, "MeV"), 1.0e-3);
  1756 + test_value(model.prefactor(), 2.0);
  1757 + test_value(model.index(), -2.1);
  1758 + test_value(model.pivot().MeV(), 100.0);
  1759 + test_value(model.lambda(), 1.0e-3);
  1760 + test_try_success();
  1761 + }
  1762 + catch (std::exception &e) {
  1763 + test_try_failure(e);
  1764 + }
  1765 +
  1766 + // Test XML constructor and value
  1767 + test_try("Test XML constructor, value and gradients");
  1768 + try {
  1769 + // Test XML constructor
  1770 + GXml xml(m_xml_model_point_eplaw2);
  1771 + GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
  1772 + GModelSpectralExpPlaw2 model(*element);
  1773 + test_value(model.size(), 4);
  1774 + test_assert(model.type() == "ExpCutoff2", "Expected \"ExpCutoff2\"");
  1775 + test_value(model.prefactor(), 5.7e-16);
  1776 + test_value(model.index(), -2.48);
  1777 + test_value(model.pivot().TeV(), 0.3);
  1778 + test_value(model.lambda(), 1.0e-6);
  1779 +
  1780 + // Test prefactor method
  1781 + model.prefactor(2.3e-16);
  1782 + test_value(model.prefactor(), 2.3e-16);
  1783 +
  1784 + // Test index method
  1785 + model.index(-2.6);
  1786 + test_value(model.index(), -2.6);
  1787 +
  1788 + // Test pivot method
  1789 + model.pivot(GEnergy(0.5, "TeV"));
  1790 + test_value(model.pivot().TeV(), 0.5);
  1791 +
  1792 + // Test cutoff parameter method
  1793 + model.lambda(4.2e-9);
  1794 + test_value(model.lambda(), 4.2e-9);
  1795 +
  1796 + // Test operator access
  1797 + const char* strarray[] = {"Prefactor", "Index", "PivotEnergy", "lambda"};
  1798 + for (int i = 0; i < 4; ++i) {
  1799 + std::string keyname(strarray[i]);
  1800 + model[keyname].remove_range(); // To allow setting of any value
  1801 + model[keyname].value(2.1);
  1802 + model[keyname].error(1.9);
  1803 + model[keyname].gradient(0.8);
  1804 + test_value(model[keyname].value(), 2.1);
  1805 + test_value(model[keyname].error(), 1.9);
  1806 + test_value(model[keyname].gradient(), 0.8);
  1807 + }
  1808 +
  1809 + // Success if we reached this point
  1810 + test_try_success();
  1811 + }
  1812 + catch (std::exception &e) {
  1813 + test_try_failure(e);
  1814 + }
  1815 +
  1816 + // Exit test
  1817 + return;
  1818 +}
  1819 +
  1820 +
  1821 +/***********************************************************************//**
1734 1822 * @brief Test GModelSpectralSuperExpPlaw class
1735 1823 ***************************************************************************/
1736 1824 void TestGModel::test_supeplaw(void)
... ...
test/test_GModel.hpp
  View file @52e5c22
... ... @@ -64,6 +64,7 @@ public:
64 64 void test_plaw(void);
65 65 void test_plaw2(void);
66 66 void test_eplaw(void);
  67 + void test_eplaw2(void);
67 68 void test_bplaw(void);
68 69 void test_supeplaw(void);
69 70 void test_logparabola(void);
... ... @@ -88,6 +89,7 @@ private:
88 89 std::string m_xml_model_point_plaw;
89 90 std::string m_xml_model_point_plaw2;
90 91 std::string m_xml_model_point_eplaw;
  92 + std::string m_xml_model_point_eplaw2;
91 93 std::string m_xml_model_point_bplaw;
92 94 std::string m_xml_model_point_supeplaw;
93 95 std::string m_xml_model_point_logparabola;
... ...