Dmitry Malyshev / gammalib

Commit 1006ac2be6b122d4d1ec15049edb8c29536bf727

Authored by ellisowen
1 parent 3a61e2e2

Add GModelSpectralGauss class (issue #1014)

Inline Side-by-side
Showing 8 changed files with 1043 additions and 1 deletions
... ... @@ -2,7 +2,7 @@ New Features and Important Changes in GammaLib 00-08-00
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -9 January 2014
  5 +13 January 2014
6 6  
7 7  
8 8 1. Introduction
... ... @@ -496,6 +496,7 @@ GOptimizerPar are the XML read() and write() methods.
496 496 A broken power law spectral model has been implemented by the
497 497 GModelSpectralBrokenPlaw class.
498 498  
  499 +A Gauss spectral model has been implemented by the GModelSpectralGauss class.
499 500  
500 501 9. Numerics module
501 502 ------------------
... ...
include/GModelSpectralGauss.hpp 0 → 100644
  View file @1006ac2
  1 +/***************************************************************************
  2 + * GModelSpectralGauss.hpp - Spectral Gaussian model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2009-2013 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 GModelSpectralGauss.hpp
  23 + * @brief Gaussian spectral model class interface definition
  24 + * @author Christoph Deil & Ellis Owen
  25 + */
  26 +
  27 +#ifndef GMODELSPECTRALGAUSS_HPP
  28 +#define GMODELSPECTRALGAUSS_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GModelPar.hpp"
  33 +#include "GModelSpectral.hpp"
  34 +#include "GEnergy.hpp"
  35 +#include "GXmlElement.hpp"
  36 +#include "GFunction.hpp"
  37 +
  38 +
  39 +/***********************************************************************//**
  40 + * @class GModelSpectralGauss
  41 + *
  42 + * @brief Gaussian spectral model class
  43 + *
  44 + * This class implements a Gaussian spectrum. The model is defined by
  45 + *
  46 + * \f[
  47 + * S_{\rm E}(E | t) =
  48 + * \frac{\tt m\_norm}{\sqrt{2 \pi} \sigma^2} \exp
  49 + * \left(\frac{-(E - m\_mean)^2}{2 m_\sigma^2} \right)
  50 + * \f]
  51 + *
  52 + * where
  53 + * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  54 + * - \f${\tt m\_mean}\f$ is the mean energy,
  55 + * - \f${\tt m\_sigma}\f$ is the energy width.
  56 + ***************************************************************************/
  57 +class GModelSpectralGauss : public GModelSpectral {
  58 +
  59 +public:
  60 + // Constructors and destructors
  61 + GModelSpectralGauss(void);
  62 + explicit GModelSpectralGauss(const GXmlElement& xml);
  63 + explicit GModelSpectralGauss(const double& prefactor,
  64 + const GEnergy& mean,
  65 + const GEnergy& sigma);
  66 + GModelSpectralGauss(const GModelSpectralGauss& model);
  67 + virtual ~GModelSpectralGauss(void);
  68 +
  69 + // Operators
  70 + virtual GModelSpectralGauss& operator=(const GModelSpectralGauss& model);
  71 +
  72 + // Implemented pure virtual base class methods
  73 + virtual void clear(void);
  74 + virtual GModelSpectralGauss* clone(void) const;
  75 + virtual std::string type(void) const;
  76 + virtual double eval(const GEnergy& srcEng,
  77 + const GTime& srcTime) const;
  78 + virtual double eval_gradients(const GEnergy& srcEng,
  79 + const GTime& srcTime);
  80 + virtual double flux(const GEnergy& emin,
  81 + const GEnergy& emax) const;
  82 + virtual double eflux(const GEnergy& emin,
  83 + const GEnergy& emax) const;
  84 + virtual GEnergy mc(const GEnergy& emin,
  85 + const GEnergy& emax,
  86 + const GTime& time,
  87 + GRan& ran) const;
  88 + virtual void read(const GXmlElement& xml);
  89 + virtual void write(GXmlElement& xml) const;
  90 + virtual std::string print(const GChatter& chatter = NORMAL) const;
  91 +
  92 + // Other methods
  93 + double prefactor(void) const;
  94 + void prefactor(const double& prefactor);
  95 + GEnergy mean(void) const;
  96 + void mean(const GEnergy& mean);
  97 + GEnergy sigma(void) const;
  98 + void sigma(const GEnergy& sigma);
  99 +
  100 +
  101 +
  102 +protected:
  103 + // Protected methods
  104 + void init_members(void);
  105 + void copy_members(const GModelSpectralGauss& model);
  106 + void free_members(void);
  107 +
  108 + // Energy flux integration kernel
  109 + class eflux_kernel : public GFunction {
  110 + public:
  111 + eflux_kernel(const double& norm,
  112 + const double& mean,
  113 + const double& sigma) :
  114 + m_norm(norm),
  115 + m_mean(mean),
  116 + m_sigma(sigma) {}
  117 + double eval(const double& eng);
  118 + protected:
  119 + const double& m_norm; //!< Normalization
  120 + const double& m_mean; //!< Mean
  121 + const double& m_sigma; //!< Sigma
  122 + };
  123 +
  124 + // Protected members
  125 + GModelPar m_norm; //!< Normalization factor
  126 + GModelPar m_mean; //!< Gaussian mean energy
  127 + GModelPar m_sigma; //!< Gaussian energy width
  128 +};
  129 +
  130 +
  131 +/***********************************************************************//**
  132 + * @brief Return model type
  133 + *
  134 + * @return "ConstantValue".
  135 + *
  136 + * Returns the type of the Gaussian spectral model.
  137 + ***************************************************************************/
  138 +inline
  139 +std::string GModelSpectralGauss::type(void) const
  140 +{
  141 + return "Gaussian";
  142 +}
  143 +
  144 +
  145 +/***********************************************************************//**
  146 + * @brief Return pre factor
  147 + *
  148 + * @return Pre factor (ph/cm2/s/MeV).
  149 + ***************************************************************************/
  150 +inline
  151 +double GModelSpectralGauss::prefactor(void) const
  152 +{
  153 + return (m_norm.value());
  154 +}
  155 +
  156 +
  157 +/***********************************************************************//**
  158 + * @brief Set pre factor
  159 + *
  160 + * @param[in] prefactor Pre factor (ph/cm2/s/MeV).
  161 + ***************************************************************************/
  162 +inline
  163 +void GModelSpectralGauss::prefactor(const double& prefactor)
  164 +{
  165 + m_norm.value(prefactor);
  166 + return;
  167 +}
  168 +
  169 +/***********************************************************************//**
  170 + * @brief Return mean energy
  171 + *
  172 + * @return Mean energy.
  173 + ***************************************************************************/
  174 +inline
  175 +GEnergy GModelSpectralGauss::mean(void) const
  176 +{
  177 + GEnergy energy;
  178 + energy.MeV(m_mean.value());
  179 + return energy;
  180 +}
  181 +
  182 +
  183 +/***********************************************************************//**
  184 + * @brief Set mean energy
  185 + *
  186 + * @param[in] Mean energy.
  187 + ***************************************************************************/
  188 +inline
  189 +void GModelSpectralGauss::mean(const GEnergy& mean)
  190 +{
  191 + m_mean.value(mean.MeV());
  192 + return;
  193 +}
  194 +
  195 +/***********************************************************************//**
  196 + * @brief Return energy width
  197 + *
  198 + * @return Energy width
  199 + ***************************************************************************/
  200 +inline
  201 +GEnergy GModelSpectralGauss::sigma(void) const
  202 +{
  203 + GEnergy energy;
  204 + energy.MeV(m_sigma.value());
  205 + return energy;
  206 +}
  207 +
  208 +
  209 +/***********************************************************************//**
  210 + * @brief Set energy width
  211 + *
  212 + * @param[in] Energy width
  213 + ***************************************************************************/
  214 +inline
  215 +void GModelSpectralGauss::sigma(const GEnergy& sigma)
  216 +{
  217 + m_sigma.value(sigma.MeV());
  218 + return;
  219 +}
  220 +
  221 +
  222 +#endif /* GMODELSPECTRALGAUSS_HPP */
... ...
include/GammaLib.hpp
  View file @1006ac2
... ... @@ -203,6 +203,7 @@
203 203 #include "GModelSpectralConst.hpp"
204 204 #include "GModelSpectralExpPlaw.hpp"
205 205 #include "GModelSpectralFunc.hpp"
  206 +#include "GModelSpectralGauss.hpp"
206 207 #include "GModelSpectralLogParabola.hpp"
207 208 #include "GModelSpectralNodes.hpp"
208 209 #include "GModelSpectralPlaw.hpp"
... ...
include/Makefile.am
  View file @1006ac2
... ... @@ -140,6 +140,7 @@ pkginclude_HEADERS = GammaLib.hpp \
140 140 GModelSpectralConst.hpp \
141 141 GModelSpectralExpPlaw.hpp \
142 142 GModelSpectralFunc.hpp \
  143 + GModelSpectralGauss.hpp \
143 144 GModelSpectralLogParabola.hpp \
144 145 GModelSpectralNodes.hpp \
145 146 GModelSpectralPlaw.hpp \
... ...
pyext/GModelSpectralGauss.i 0 → 100644
  View file @1006ac2
  1 +/***************************************************************************
  2 + * GModelSpectralGauss.i - Spectral gaussian model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2011-2013 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 GModelSpectralGauss.i
  23 + * @brief Spectral constant model class interface definition
  24 + * @author Ellis Owen, Christoph Deil
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelSpectralGauss.hpp"
  29 +#include "GTools.hpp"
  30 +#include "GEnergy.hpp"
  31 +%}
  32 +
  33 +
  34 +/***********************************************************************//**
  35 + * @class GModelSpectralGauss
  36 + *
  37 + * @brief Spectral gaussian model class
  38 + ***************************************************************************/
  39 +class GModelSpectralGauss : public GModelSpectral {
  40 +public:
  41 + // Constructors and destructors
  42 + GModelSpectralGauss(void);
  43 + explicit GModelSpectralGauss(const GXmlElement& xml);
  44 + explicit GModelSpectralGauss(const double& prefactor,
  45 + const GEnergy& mean,
  46 + const GEnergy& sigma);
  47 + GModelSpectralGauss(const GModelSpectralGauss& model);
  48 + virtual ~GModelSpectralGauss(void);
  49 +
  50 + // Implemented pure virtual base class methods
  51 + virtual void clear(void);
  52 + virtual GModelSpectralGauss* clone(void) const;
  53 + virtual std::string type(void) const;
  54 + virtual double eval(const GEnergy& srcEng,
  55 + const GTime& srcTime) const;
  56 + virtual double eval_gradients(const GEnergy& srcEng,
  57 + const GTime& srcTime);
  58 + virtual double flux(const GEnergy& emin,
  59 + const GEnergy& emax) const;
  60 + virtual double eflux(const GEnergy& emin,
  61 +
  62 + const GEnergy& emax) const;
  63 + virtual GEnergy mc(const GEnergy& emin,
  64 + const GEnergy& emax,
  65 + const GTime& time,
  66 + GRan& ran) const;
  67 + virtual void read(const GXmlElement& xml);
  68 + virtual void write(GXmlElement& xml) const;
  69 +
  70 + // Other methods
  71 + double prefactor(void) const;
  72 + void prefactor(const double& prefactor);
  73 + GEnergy mean(void) const;
  74 + void mean(const GEnergy& mean);
  75 + GEnergy sigma(void) const;
  76 + void sigma(const GEnergy& sigma);
  77 +};
  78 +
  79 +
  80 +/***********************************************************************//**
  81 + * @brief GModelSpectralGauss class extension
  82 + ***************************************************************************/
  83 +%extend GModelSpectralGauss {
  84 + GModelSpectralGauss copy() {
  85 + return (*self);
  86 + }
  87 +};
... ...
pyext/gammalib/model.i
  View file @1006ac2
... ... @@ -77,6 +77,7 @@
77 77 %include "GModelSpectralConst.i"
78 78 %include "GModelSpectralExpPlaw.i"
79 79 %include "GModelSpectralFunc.i"
  80 +%include "GModelSpectralGauss.i"
80 81 %include "GModelSpectralLogParabola.i"
81 82 %include "GModelSpectralNodes.i"
82 83 %include "GModelSpectralPlaw.i"
... ...
src/model/GModelSpectralGauss.cpp 0 → 100644
  View file @1006ac2
  1 +/***************************************************************************
  2 + * GModelSpectralGauss.cpp - Spectral Gaussian model class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2009-2013 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 GModelSpectralGauss.cpp
  23 + * @brief Gaussian spectral model class implementation
  24 + * @author Christoph Deil & Ellis Owen
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "GException.hpp"
  32 +#include "GTools.hpp"
  33 +#include "GMath.hpp"
  34 +#include "GModelSpectralGauss.hpp"
  35 +#include "GModelSpectralRegistry.hpp"
  36 +#include "GFunction.hpp"
  37 +#include "GIntegral.hpp"
  38 +
  39 +/* __ Constants __________________________________________________________ */
  40 +
  41 +/* __ Globals ____________________________________________________________ */
  42 +const GModelSpectralGauss g_spectral_gauss_seed;
  43 +const GModelSpectralRegistry g_spectral_gauss_registry(&g_spectral_gauss_seed);
  44 +
  45 +/* __ Method name definitions ____________________________________________ */
  46 +#define G_FLUX "GModelSpectralGauss::flux(GEnergy&, GEnergy&)"
  47 +#define G_EFLUX "GModelSpectralGauss::eflux(GEnergy&, GEnergy&)"
  48 +#define G_MC "GModelSpectralGauss::mc(GEnergy&, GEnergy&, GRan&)"
  49 +#define G_READ "GModelSpectralGauss::read(GXmlElement&)"
  50 +#define G_WRITE "GModelSpectralGauss::write(GXmlElement&)"
  51 +
  52 +/* __ Macros _____________________________________________________________ */
  53 +
  54 +/* __ Coding definitions _________________________________________________ */
  55 +
  56 +/* __ Debug definitions __________________________________________________ */
  57 +
  58 +
  59 +/*==========================================================================
  60 + = =
  61 + = Constructors/destructors =
  62 + = =
  63 + ==========================================================================*/
  64 +
  65 +/***********************************************************************//**
  66 + * @brief Void constructor
  67 + ***************************************************************************/
  68 +GModelSpectralGauss::GModelSpectralGauss(void) : GModelSpectral()
  69 +{
  70 + // Initialise members
  71 + init_members();
  72 +
  73 + // Return
  74 + return;
  75 +}
  76 +
  77 +
  78 +/***********************************************************************//**
  79 + * @brief XML constructor
  80 + *
  81 + * @param[in] xml XML element.
  82 + *
  83 + * Constructs constant spectral model by extracting information from an XML
  84 + * element. See the read() method for more information about the expected
  85 + * structure of the XML element.
  86 + ***************************************************************************/
  87 +GModelSpectralGauss::GModelSpectralGauss(const GXmlElement& xml) :
  88 + GModelSpectral()
  89 +{
  90 + // Initialise members
  91 + init_members();
  92 +
  93 + // Read information from XML element
  94 + read(xml);
  95 +
  96 + // Return
  97 + return;
  98 +}
  99 +
  100 +
  101 +/***********************************************************************//**
  102 + * @brief Constructor
  103 + *
  104 + * @param[in] prefactor Power law pre factor (in ph/cm2/s/MeV).
  105 + * @param[in] mean Mean energy.
  106 + * @param[in] sigma Energy width.
  107 + ***************************************************************************/
  108 +GModelSpectralGauss::GModelSpectralGauss(const double& prefactor,
  109 + const GEnergy& mean,
  110 + const GEnergy& sigma) :
  111 + GModelSpectral()
  112 +{
  113 + // Initialise members
  114 + init_members();
  115 +
  116 + // Set parameters
  117 + m_norm.value(prefactor);
  118 + m_mean.value(mean.MeV());
  119 + m_sigma.value(sigma.MeV());
  120 +
  121 + // Autoscale parameters
  122 + autoscale();
  123 +
  124 + // Return
  125 + return;
  126 +}
  127 +
  128 +
  129 +/***********************************************************************//**
  130 + * @brief Copy constructor
  131 + *
  132 + * @param[in] model Spectral Gaussian model.
  133 + ***************************************************************************/
  134 +GModelSpectralGauss::GModelSpectralGauss(const GModelSpectralGauss& model) :
  135 + GModelSpectral(model)
  136 +{
  137 + // Initialise members
  138 + init_members();
  139 +
  140 + // Copy members
  141 + copy_members(model);
  142 +
  143 + // Return
  144 + return;
  145 +}
  146 +
  147 +
  148 +/***********************************************************************//**
  149 + * @brief Destructor
  150 + ***************************************************************************/
  151 +GModelSpectralGauss::~GModelSpectralGauss(void)
  152 +{
  153 + // Free members
  154 + free_members();
  155 +
  156 + // Return
  157 + return;
  158 +}
  159 +
  160 +
  161 +/*==========================================================================
  162 + = =
  163 + = Operators =
  164 + = =
  165 + ==========================================================================*/
  166 +
  167 +/***********************************************************************//**
  168 + * @brief Assignment operator
  169 + *
  170 + * @param[in] model Gauss spectral model.
  171 + * @return Gauss spectral model.
  172 + ***************************************************************************/
  173 +GModelSpectralGauss& GModelSpectralGauss::operator=(const GModelSpectralGauss& model)
  174 +{
  175 + // Execute only if object is not identical
  176 + if (this != &model) {
  177 +
  178 + // Copy base class members
  179 + this->GModelSpectral::operator=(model);
  180 +
  181 + // Free members
  182 + free_members();
  183 +
  184 + // Initialise members
  185 + init_members();
  186 +
  187 + // Copy members
  188 + copy_members(model);
  189 +
  190 + } // endif: object was not identical
  191 +
  192 + // Return
  193 + return *this;
  194 +}
  195 +
  196 +
  197 +/*==========================================================================
  198 + = =
  199 + = Public methods =
  200 + = =
  201 + ==========================================================================*/
  202 +
  203 +/***********************************************************************//**
  204 + * @brief Clear Gaussian spectral model
  205 + ***************************************************************************/
  206 +void GModelSpectralGauss::clear(void)
  207 +{
  208 + // Free class members (base and derived classes, derived class first)
  209 + free_members();
  210 + this->GModelSpectral::free_members();
  211 +
  212 + // Initialise members
  213 + this->GModelSpectral::init_members();
  214 + init_members();
  215 +
  216 + // Return
  217 + return;
  218 +}
  219 +
  220 +
  221 +/***********************************************************************//**
  222 + * @brief Clone Gaussian spectral model
  223 + *
  224 + * @return Pointer to deep copy of Gaussian spectral model.
  225 + ***************************************************************************/
  226 +GModelSpectralGauss* GModelSpectralGauss::clone(void) const
  227 +{
  228 + // Clone Gaussian spectral model
  229 + return new GModelSpectralGauss(*this);
  230 +}
  231 +
  232 +
  233 +/***********************************************************************//**
  234 + * @brief Evaluate model value
  235 + *
  236 + * @param[in] srcEng True photon energy.
  237 + * @param[in] srcTime True photon arrival time.
  238 + * @return Model value (ph/cm2/s/MeV).
  239 + *
  240 + * Evaluates
  241 + *
  242 + * \f[
  243 + * \frac{dN}{dE}=\frac{\norm}{\sqrt{2\pi}\sigma}\exp(\frac{(E-\bar{E})^2}{2\sigma^2})
  244 + * \f]
  245 + ***************************************************************************/
  246 +double GModelSpectralGauss::eval(const GEnergy& srcEng,
  247 + const GTime& srcTime) const
  248 +{
  249 + double energy = srcEng.MeV();
  250 + double norm = m_norm.value();
  251 + double mean = m_mean.value();
  252 + double sigma = m_sigma.value();
  253 +
  254 + // Compute function value
  255 + double term1 = (norm / sigma) * gammalib::inv_sqrt2pi;
  256 + double term2 = (energy - mean) * (energy - mean) / (2 * sigma * sigma);
  257 + double value = term1 * std::exp(- term2);
  258 +
  259 + // Return
  260 + return value;
  261 +}
  262 +
  263 +
  264 +/***********************************************************************//**
  265 + * @brief Evaluate model value and gradient
  266 + *
  267 + * @param[in] srcEng True photon energy.
  268 + * @param[in] srcTime True photon arrival time.
  269 + * @return Model value (ph/cm2/s/MeV).
  270 + *
  271 + * This method simply calls the eval() method as no analytical gradients will
  272 + * be computed. See the eval() method for details.
  273 + ***************************************************************************/
  274 +double GModelSpectralGauss::eval_gradients(const GEnergy& srcEng,
  275 + const GTime& srcTime)
  276 +{
  277 + // Return
  278 + return eval(srcEng, srcTime);
  279 +}
  280 +
  281 +
  282 +/***********************************************************************//**
  283 + * @brief Returns model photon flux between [emin, emax] (ph/cm2/s)
  284 + *
  285 + * @param[in] emin Minimum photon energy.
  286 + * @param[in] emax Maximum photon energy.
  287 + * @return Photon flux (ph/cm2/s).
  288 + *
  289 + * Computes
  290 + *
  291 + * \f[
  292 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
  293 + * \f]
  294 + *
  295 + * where
  296 + * - [@p emin, @p emax] is an energy interval, and
  297 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  298 + * The integration is done analytically.
  299 + ***************************************************************************/
  300 +
  301 +double GModelSpectralGauss::flux(const GEnergy& emin,
  302 + const GEnergy& emax) const{
  303 +
  304 + double energy_min = emin.MeV();
  305 + double energy_max = emax.MeV();
  306 +
  307 +
  308 + double norm = m_norm.value();
  309 + double mean = m_mean.value();
  310 + double sigma = m_sigma.value();
  311 +
  312 + double zmin = (energy_min - mean)/sqrt(2.)/sigma;
  313 + double zmax = (energy_max - mean)/sqrt(2.)/sigma;
  314 +
  315 + double flux = 0.0;
  316 +
  317 + // Compute only if integration range is valid
  318 + if (emin < emax) {
  319 +
  320 + // Compute flux for a constant model
  321 + flux = norm*(gammalib::erfcc(zmin) - gammalib::erfcc(zmax))/2.;
  322 +
  323 + } // endif: integration range was valid
  324 +
  325 + // Return
  326 + return flux;
  327 +}
  328 +
  329 +
  330 +/***********************************************************************//**
  331 + * @brief Returns model energy flux between [emin, emax] (erg/cm2/s)
  332 + *
  333 + * @param[in] emin Minimum photon energy.
  334 + * @param[in] emax Maximum photon energy.
  335 + * @return Energy flux (erg/cm2/s).
  336 + *
  337 + * Computes
  338 + *
  339 + * \f[
  340 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
  341 + * \f]
  342 + *
  343 + * where
  344 + * - [@p emin, @p emax] is an energy interval, and
  345 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  346 + * The integration is done analytically.
  347 + ***************************************************************************/
  348 +double GModelSpectralGauss::eflux(const GEnergy& emin,
  349 + const GEnergy& emax) const
  350 +{
  351 + // Initialise flux
  352 + double eflux = 0.0;
  353 +
  354 + // Compute only if integration range is valid
  355 + if (emin < emax) {
  356 +
  357 + // Setup integration kernel
  358 + eflux_kernel integrand(m_norm.value(), m_mean.value(),
  359 + m_sigma.value());
  360 + GIntegral integral(&integrand);
  361 +
  362 + // Get integration boundaries in MeV
  363 + double e_min = emin.MeV();
  364 + double e_max = emax.MeV();
  365 +
  366 + // Perform integration
  367 + eflux = integral.romb(e_min, e_max);
  368 +
  369 + // Convert from MeV/cm2/s to erg/cm2/s
  370 + eflux *= gammalib::MeV2erg;
  371 +
  372 + } // endif: integration range was valid
  373 +
  374 + // Return
  375 + return eflux;
  376 +}
  377 +
  378 +
  379 +/***********************************************************************//**
  380 + * @brief Returns MC energy between [emin, emax]
  381 + *
  382 + * @param[in] emin Minimum photon energy.
  383 + * @param[in] emax Maximum photon energy.
  384 + * @param[in] time True photon arrival time.
  385 + * @param[in,out] ran Random number generator.
  386 + * @return Energy.
  387 + *
  388 + * @exception GException::erange_invalid
  389 + * Energy range is invalid (emin < emax required).
  390 + *
  391 + * Returns Monte Carlo energy by randomly drawing from a constant between
  392 + * the minimum and maximum photon energy.
  393 + *
  394 + * Method Used: Box-Muller transform, outlined here:
  395 + * http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
  396 + *
  397 + * Code from: http://www.design.caltech.edu/erik/Misc/Gaussian.html
  398 + ***************************************************************************/
  399 +GEnergy GModelSpectralGauss::mc(const GEnergy& emin,
  400 + const GEnergy& emax,
  401 + const GTime& time,
  402 + GRan& ran) const
  403 +{
  404 + double xmax = emax.MeV();
  405 + double xmin = emin.MeV();
  406 + double energy;
  407 +
  408 + // Throw an exception if energy range is invalid
  409 + if (xmin >= xmax) {
  410 + throw GException::erange_invalid(G_MC, xmin, xmax,
  411 + "Minimum energy < maximum energy required.");
  412 + }
  413 +
  414 + // Get uniform value between 0 and 1
  415 + double u = ran.uniform();
  416 + float x1, x2, w;
  417 +
  418 + do {
  419 + do{
  420 + x1 = 2.0 * u - 1.0;
  421 + x2 = 2.0 * u - 1.0;
  422 + w = x1 * x1 + x2 * x2;
  423 + } while (w >= 1.0);
  424 +
  425 + w = sqrt( (-2.0 * log( w ) ) / w );
  426 +
  427 + double val = x1 * w;
  428 +
  429 + // Map into [emin,emax] range
  430 +
  431 + energy = (xmax - xmin) * val + xmin;
  432 +
  433 + } while ((xmin <= energy) && (energy < xmax));
  434 +
  435 + // Return energy
  436 + return GEnergy(energy, "MeV");
  437 +}
  438 +
  439 +
  440 +/***********************************************************************//**
  441 + * @brief Read model from XML element
  442 + *
  443 + * @param[in] xml XML element containing Gaussian model information.
  444 + *
  445 + * @exception GException::model_invalid_parnum
  446 + * Invalid number of model parameters found in XML element.
  447 + * @exception GException::model_invalid_parnames
  448 + * Invalid model parameter name found in XML element.
  449 + *
  450 + * Read the spectral Gaussian information from an XML element.
  451 + * The format of the XML elements is:
  452 + *
  453 + * <spectrum type="Gaussian">
  454 + * <parameter name="Prefactor" scale=".." value=".." min=".." max=".." free=".."/>
  455 + * <parameter name="Mean" scale=".." value=".." min=".." max=".." free=".."/>
  456 + * <parameter name="Sigma" scale=".." value=".." min=".." max=".." free=".."/>
  457 + * </spectrum>
  458 + *
  459 + ***************************************************************************/
  460 +void GModelSpectralGauss::read(const GXmlElement& xml)
  461 +{
  462 + const int n_pars = 3;
  463 +
  464 + // Verify that XML element has exactly 3 parameters
  465 + if (xml.elements() != n_pars || xml.elements("parameter") != n_pars) {
  466 + throw GException::model_invalid_parnum(G_READ, xml,
  467 + "Gaussian model requires exactly 3 parameters.");
  468 + }
  469 +
  470 + // Extract model parameters
  471 + int npar[] = {0, 0, 0};
  472 + for (int i = 0; i < n_pars; ++i) {
  473 +
  474 + // Get parameter element
  475 + const GXmlElement* par = xml.element("parameter", i);
  476 +
  477 + // Handle prefactor
  478 + if (par->attribute("name") == "Prefactor") {
  479 + m_norm.read(*par);
  480 + npar[0]++;
  481 + }
  482 +
  483 + // Handle mean
  484 + else if (par->attribute("name") == "Mean") {
  485 + m_mean.read(*par);
  486 + npar[1]++;
  487 + }
  488 +
  489 + // Handle sigma
  490 + else if (par->attribute("name") == "Sigma") {
  491 + m_sigma.read(*par);
  492 + npar[2]++;
  493 + }
  494 +
  495 + } // endfor: looped over all parameters
  496 +
  497 + // Verify that all parameters were found
  498 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1) {
  499 + throw GException::model_invalid_parnames(G_READ, xml,
  500 + "Require \"Prefactor\", \"Mean\", and \"Sigma\""
  501 + " parameters.");
  502 + }
  503 +
  504 + // Return
  505 + return;
  506 +}
  507 +
  508 +
  509 +/***********************************************************************//**
  510 + * @brief Write model into XML element
  511 + *
  512 + * @param[in] xml XML element into which model information is written.
  513 + *
  514 + * @exception GException::model_invalid_spectral
  515 + * Existing XML element is not of type "Gaussian"
  516 + * @exception GException::model_invalid_parnum
  517 + * Invalid number of model parameters or nodes found in XML element.
  518 + * @exception GException::model_invalid_parnames
  519 + * Invalid model parameter names found in XML element.
  520 + *
  521 + * Writes the spectral information into an XML element. The format of the XML
  522 + * element is
  523 + *
  524 + * <spectrum type="Gaussian">
  525 + * <parameter name="Prefactor" scale=".." value=".." min=".." max=".." free=".."/>
  526 + * <parameter name="Mean" scale=".." value=".." min=".." max=".." free=".."/>
  527 + * <parameter name="Sigma" scale=".." value=".." min=".." max=".." free=".."/>
  528 + * </spectrum>
  529 + ***************************************************************************/
  530 +void GModelSpectralGauss::write(GXmlElement& xml) const
  531 +{
  532 + // Set model type
  533 + if (xml.attribute("type") == "") {
  534 + xml.attribute("type", "Gaussian");
  535 + }
  536 +
  537 + // Verify model type
  538 + if (xml.attribute("type") != "Gaussian") {
  539 + throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"),
  540 + "Spectral model is not of type \"Gaussian\".");
  541 + }
  542 +
  543 + // If XML element has 0 nodes then append 3 parameter nodes
  544 + if (xml.elements() == 0) {
  545 + xml.append(GXmlElement("parameter name=\"Prefactor\""));
  546 + xml.append(GXmlElement("parameter name=\"Mean\""));
  547 + xml.append(GXmlElement("parameter name=\"Sigma\""));
  548 + }
  549 +
  550 + // Verify that XML element has exactly 3 parameters
  551 + if (xml.elements() != 3 || xml.elements("parameter") != 3) {
  552 + throw GException::model_invalid_parnum(G_WRITE, xml,
  553 + "Power law model requires exactly 3 parameters.");
  554 + }
  555 +
  556 + // Set or update model parameter attributes
  557 + int npar[] = {0, 0, 0};
  558 + for (int i = 0; i < 3; ++i) {
  559 +
  560 + // Get parameter element
  561 + GXmlElement* par = xml.element("parameter", i);
  562 +
  563 + // Handle prefactor
  564 + if (par->attribute("name") == "Prefactor") {
  565 + npar[0]++;
  566 + m_norm.write(*par);
  567 + }
  568 +
  569 + // Handle mean
  570 + else if (par->attribute("name") == "Mean") {
  571 + npar[1]++;
  572 + m_mean.write(*par);
  573 + }
  574 +
  575 + // Handle sigma
  576 + else if (par->attribute("name") == "Sigma") {
  577 + npar[2]++;
  578 + m_sigma.write(*par);
  579 + }
  580 +
  581 +
  582 + } // endfor: looped over all parameters
  583 +
  584 + // Check of all required parameters are present
  585 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 ) {
  586 + throw GException::model_invalid_parnames(G_WRITE, xml,
  587 + "Require \"Prefactor\", \"Mean\" and \"Sigma\""
  588 + " parameters.");
  589 + }
  590 +
  591 + // Return
  592 + return;
  593 +}
  594 +
  595 +
  596 +/***********************************************************************//**
  597 + * @brief Print spectral model information
  598 + *
  599 + * @param[in] chatter Chattiness (defaults to NORMAL).
  600 + * @return String containing spectral model information.
  601 + ***************************************************************************/
  602 +std::string GModelSpectralGauss::print(const GChatter& chatter) const
  603 +{
  604 + // Initialise result string
  605 + std::string result;
  606 +
  607 + // Continue only if chatter is not silent
  608 + if (chatter != SILENT) {
  609 +
  610 + // Append header
  611 + result.append("=== GModelSpectralGauss ===");
  612 +
  613 + // Append model content
  614 + result.append("\n"+gammalib::parformat("Number of parameters"));
  615 + result.append(gammalib::str(size()));
  616 + for (int i = 0; i < size(); ++i) {
  617 + result.append("\n"+m_pars[i]->print(chatter));
  618 + }
  619 +
  620 + } // endif: chatter was not silent
  621 +
  622 + // Return result
  623 + return result;
  624 +}
  625 +
  626 +
  627 +/*==========================================================================
  628 + = =
  629 + = Private methods =
  630 + = =
  631 + ==========================================================================*/
  632 +
  633 +/***********************************************************************//**
  634 + * @brief Initialise class members
  635 + ***************************************************************************/
  636 +void GModelSpectralGauss::init_members(void)
  637 +{
  638 + // Initialise normalisation
  639 + m_norm.clear();
  640 + m_norm.name("Prefactor");
  641 + m_norm.unit("ph/cm2/s/MeV");
  642 + m_norm.scale(1.0);
  643 + m_norm.value(1.0); // default: 1.0
  644 + m_norm.min(0.0); // min: 0.0
  645 + m_norm.free();
  646 + m_norm.gradient(0.0);
  647 + m_norm.has_grad(false);
  648 +
  649 + // Initialise mean energy
  650 + m_mean.clear();
  651 + m_mean.name("Mean");
  652 + m_mean.unit("MeV");
  653 + m_mean.scale(1.0);
  654 + m_mean.value(1000.0); // default: 1000.0
  655 + m_mean.min(0.1); // min: 0.1
  656 + m_mean.free();
  657 + m_mean.gradient(0.0);
  658 + m_mean.has_grad(false);
  659 +
  660 + // Initialise energy width
  661 + m_sigma.clear();
  662 + m_sigma.name("Sigma");
  663 + m_sigma.unit("MeV");
  664 + m_sigma.scale(1.0);
  665 + m_sigma.value(1000.0); // default: 1000.0
  666 + m_sigma.min(0.1); // min: 0.1
  667 + m_sigma.free();
  668 + m_sigma.gradient(0.0);
  669 + m_sigma.has_grad(false);
  670 +
  671 + // Set parameter pointer(s)
  672 + m_pars.clear();
  673 + m_pars.push_back(&m_norm);
  674 + m_pars.push_back(&m_mean);
  675 + m_pars.push_back(&m_sigma);
  676 +
  677 + // Return
  678 + return;
  679 +}
  680 +
  681 +
  682 +/***********************************************************************//**
  683 + * @brief Copy class members
  684 + *
  685 + * @param[in] model Spectral constant model.
  686 + ***************************************************************************/
  687 +void GModelSpectralGauss::copy_members(const GModelSpectralGauss& model)
  688 +{
  689 + // Copy members
  690 + m_norm = model.m_norm;
  691 + m_mean = model.m_mean;
  692 + m_sigma = model.m_sigma;
  693 +
  694 + // Set parameter pointer(s)
  695 + m_pars.clear();
  696 + m_pars.push_back(&m_norm);
  697 + m_pars.push_back(&m_mean);
  698 + m_pars.push_back(&m_sigma);
  699 +
  700 + // Return
  701 + return;
  702 +}
  703 +
  704 +
  705 +/***********************************************************************//**
  706 + * @brief Delete class members
  707 + ***************************************************************************/
  708 +void GModelSpectralGauss::free_members(void)
  709 +{
  710 + // Return
  711 + return;
  712 +}
  713 +
  714 +/***********************************************************************//**
  715 + * @brief Kernel for energy flux integration
  716 + *
  717 + * @param[in] energy Energy (MeV).
  718 + ***************************************************************************/
  719 +double GModelSpectralGauss::eflux_kernel::eval(const double& energy)
  720 +{
  721 + // Evaluate function value
  722 + double term1 = (m_norm / m_sigma) * gammalib::inv_sqrt2pi;
  723 + double term2 = (energy - m_mean) * (energy - m_mean) / (2 * m_sigma * m_sigma);
  724 + double value = term1 * std::exp(- term2);
  725 +
  726 + // Return value
  727 + return value;
  728 +}
... ...
src/model/Makefile.am
  View file @1006ac2
... ... @@ -31,6 +31,7 @@ sources = GModelPar.cpp \
31 31 GModelSpectralConst.cpp \
32 32 GModelSpectralExpPlaw.cpp \
33 33 GModelSpectralFunc.cpp \
  34 + GModelSpectralGauss.cpp \
34 35 GModelSpectralLogParabola.cpp \
35 36 GModelSpectralNodes.cpp \
36 37 GModelSpectralPlaw.cpp \
... ...