Simon Bonnefoy / gammalib

Commit 1e21288db570b03824179bbf22035e077f068ae3

Authored by Jürgen Knödlseder
2 parents fbd66c42 915a6d61

Integrate GModelSpectralExpInvPlaw class

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