Rolf Buehler / gammalib

Commit 49606587aa03338cb45cf49abc05d9ec40db4159

Authored by Jürgen Knödlseder
2 parents ab11733f ed7891e8

Merge branch '1009-add_broken_plaw' into integration 

Conflicts:
	doc/source/user_manual/modules/model.rst
Inline Side-by-side
Showing 11 changed files with 1554 additions and 9 deletions
include/GModelSpectralBrokenPlaw.hpp 0 → 100644
  View file @4960658
  1 +/***************************************************************************
  2 + * GModelSpectralBrokenPlaw.hpp - Broken power law spectrum class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Anneli Schulz *
  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 GModelSpectralBrokenPlaw.hpp
  23 + * @brief Broken power law spectrum class definition
  24 + * @author Anneli Schulz
  25 + */
  26 +
  27 +#ifndef GMODELSPECTRALBROKENPLAW_HPP
  28 +#define GMODELSPECTRALBROKENPLAW_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <string>
  32 +#include "GModelPar.hpp"
  33 +#include "GModelSpectral.hpp"
  34 +#include "GEnergy.hpp"
  35 +#include "GXmlElement.hpp"
  36 +
  37 +
  38 +/***********************************************************************//**
  39 + * @class GModelSpectralBrokenPlaw
  40 + *
  41 + * @brief Broken power law spectral model class
  42 + *
  43 + * This class implements a broken power law spectrum. The model is defined
  44 + * by
  45 + *
  46 + * \f[
  47 + * S_{\rm E}(E | t) = {\tt m\_norm}
  48 + * \left( \frac{E}{\tt m\_breakenergy} \right)^{\tt m\_index1}
  49 + * \f]
  50 + *
  51 + * where
  52 + * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  53 + * - \f${\tt m\_index1}\f$ is the spectral index1, and
  54 + * - \f${\tt m\_breakenergy}\f$ is the breakenergy energy.
  55 + ***************************************************************************/
  56 +class GModelSpectralBrokenPlaw : public GModelSpectral {
  57 +
  58 +public:
  59 + // Constructors and destructors
  60 + GModelSpectralBrokenPlaw(void);
  61 + GModelSpectralBrokenPlaw(const double& prefactor,
  62 + const double& index1,
  63 + const GEnergy& breakenergy,
  64 + const double& index2);
  65 + explicit GModelSpectralBrokenPlaw(const GXmlElement& xml);
  66 + GModelSpectralBrokenPlaw(const GModelSpectralBrokenPlaw& model);
  67 + virtual ~GModelSpectralBrokenPlaw(void);
  68 +
  69 + // Operators
  70 + virtual GModelSpectralBrokenPlaw& operator=(const GModelSpectralBrokenPlaw& model);
  71 +
  72 + // Implemented pure virtual methods
  73 + virtual void clear(void);
  74 + virtual GModelSpectralBrokenPlaw* 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 + double index1(void) const;
  95 + double index2(void) const;
  96 + GEnergy breakenergy(void) const;
  97 + void prefactor(const double& prefactor);
  98 + void index1(const double& index1);
  99 + void index2(const double& index2);
  100 + void breakenergy(const GEnergy& breakenergy);
  101 +
  102 +protected:
  103 + // Protected methods
  104 + void init_members(void);
  105 + void copy_members(const GModelSpectralBrokenPlaw& model);
  106 + void free_members(void);
  107 + void update_eval_cache(const GEnergy& energy) const;
  108 + void update_mc_cache(const GEnergy& emin, const GEnergy& emax) const;
  109 +
  110 + // Protected members
  111 + GModelPar m_norm; //!< Normalization factor
  112 + GModelPar m_index1; //!< Spectral index1
  113 + GModelPar m_index2; //!< Spectral index2
  114 + GModelPar m_breakenergy; //!< Energy of spectral break
  115 +
  116 + // Cached members used for pre-computations
  117 + mutable GEnergy m_last_energy; //!< Last energy value
  118 + mutable double m_last_index1; //!< Last index1 parameter
  119 + mutable double m_last_index2; //!< Last index1 parameter
  120 + mutable double m_last_breakenergy; //!< Last breakenergy parameter
  121 + mutable double m_last_e_norm; //!< Last E/Ebreakenergy value
  122 + mutable double m_last_log_e_norm; //!< Last ln(E/Ebreakenergy) value
  123 + mutable double m_last_power; //!< Last power value
  124 + mutable double m_mc_emin; //!< Minimum energy
  125 + mutable double m_mc_emax; //!< Maximum energy
  126 + mutable double m_mc_exponent1; //!< Exponent (index1+1)
  127 + mutable double m_mc_exponent2; //!< Exponent (index2+1)
  128 + mutable double m_mc_pow_emin; //!< Power of minimum energy
  129 + mutable double m_mc_pow_ewidth; //!< Power of energy width
  130 + mutable std::vector<double> m_mc_cum; //!< Cumulative distribution
  131 + mutable std::vector<double> m_mc_min; //!< Lower boundary for MC
  132 + mutable std::vector<double> m_mc_max; //!< Upper boundary for MC
  133 + mutable std::vector<double> m_mc_exp; //!< Exponent for MC
  134 +};
  135 +
  136 +
  137 +/***********************************************************************//**
  138 + * @brief Return model type
  139 + *
  140 + * @return "PowerLaw".
  141 + *
  142 + * Returns the type of the spectral broken power law model.
  143 + ***************************************************************************/
  144 +inline
  145 +std::string GModelSpectralBrokenPlaw::type(void) const
  146 +{
  147 + return "BrokenPowerLaw";
  148 +}
  149 +
  150 +
  151 +/***********************************************************************//**
  152 + * @brief Return pre factor
  153 + *
  154 + * @return Pre factor (ph/cm2/s/MeV).
  155 + *
  156 + * Returns the pre factor.
  157 + ***************************************************************************/
  158 +inline
  159 +double GModelSpectralBrokenPlaw::prefactor(void) const
  160 +{
  161 + return (m_norm.value());
  162 +}
  163 +
  164 +
  165 +/***********************************************************************//**
  166 + * @brief Set pre factor
  167 + *
  168 + * @param[in] prefactor Pre factor (ph/cm2/s/MeV).
  169 + *
  170 + * Sets the pre factor.
  171 + ***************************************************************************/
  172 +inline
  173 +void GModelSpectralBrokenPlaw::prefactor(const double& prefactor)
  174 +{
  175 + m_norm.value(prefactor);
  176 + return;
  177 +}
  178 +
  179 +
  180 +/***********************************************************************//**
  181 + * @brief Return power law index1
  182 + *
  183 + * @return Power law index1.
  184 + *
  185 + * Returns the power law index1.
  186 + ***************************************************************************/
  187 +inline
  188 +double GModelSpectralBrokenPlaw::index1(void) const
  189 +{
  190 + return (m_index1.value());
  191 +}
  192 +
  193 +
  194 +/***********************************************************************//**
  195 + * @brief Set power law index1
  196 + *
  197 + * @param[in] index1 Power law index1.
  198 + *
  199 + * Sets the power law index1.
  200 + ***************************************************************************/
  201 +inline
  202 +void GModelSpectralBrokenPlaw::index1(const double& index1)
  203 +{
  204 + m_index1.value(index1);
  205 + return;
  206 +}
  207 +
  208 +
  209 +/***********************************************************************//**
  210 + * @brief Return power law index2
  211 + *
  212 + * @return Power law index2.
  213 + *
  214 + * Returns the power law index2.
  215 + ***************************************************************************/
  216 +inline
  217 +double GModelSpectralBrokenPlaw::index2(void) const
  218 +{
  219 + return (m_index2.value());
  220 +}
  221 +
  222 +
  223 +/***********************************************************************//**
  224 + * @brief Set power law index2
  225 + *
  226 + * @param[in] index2 Power law index2.
  227 + *
  228 + * Sets the power law index2.
  229 + ***************************************************************************/
  230 +inline
  231 +void GModelSpectralBrokenPlaw::index2(const double& index2)
  232 +{
  233 + m_index2.value(index2);
  234 + return;
  235 +}
  236 +
  237 +
  238 +/***********************************************************************//**
  239 + * @brief Return breakenergy energy
  240 + *
  241 + * @return breakenergy energy.
  242 + *
  243 + * Returns the breakenergy energy.
  244 + ***************************************************************************/
  245 +inline
  246 +GEnergy GModelSpectralBrokenPlaw::breakenergy(void) const
  247 +{
  248 + GEnergy energy;
  249 + energy.MeV(m_breakenergy.value());
  250 + return energy;
  251 +}
  252 +
  253 +
  254 +/***********************************************************************//**
  255 + * @brief Set breakenergy energy
  256 + *
  257 + * @param[in] breakenergy breakenergy energy.
  258 + *
  259 + * Sets the breakenergy energy.
  260 + ***************************************************************************/
  261 +inline
  262 +void GModelSpectralBrokenPlaw::breakenergy(const GEnergy& breakenergy)
  263 +{
  264 + m_breakenergy.value(breakenergy.MeV());
  265 + return;
  266 +}
  267 +
  268 +#endif /* GMODELSPECTRALPLAW_HPP */
... ...
include/GammaLib.hpp
  View file @4960658
... ... @@ -199,13 +199,14 @@
199 199 #include "GModelSpatialDiffuseMap.hpp"
200 200 #include "GModelSpectral.hpp"
201 201 #include "GModelSpectralRegistry.hpp"
  202 +#include "GModelSpectralBrokenPlaw.hpp"
202 203 #include "GModelSpectralConst.hpp"
203 204 #include "GModelSpectralExpPlaw.hpp"
204 205 #include "GModelSpectralFunc.hpp"
  206 +#include "GModelSpectralLogParabola.hpp"
205 207 #include "GModelSpectralNodes.hpp"
206 208 #include "GModelSpectralPlaw.hpp"
207 209 #include "GModelSpectralPlaw2.hpp"
208   -#include "GModelSpectralLogParabola.hpp"
209 210 #include "GModelTemporal.hpp"
210 211 #include "GModelTemporalRegistry.hpp"
211 212 #include "GModelTemporalConst.hpp"
... ...
include/Makefile.am
  View file @4960658
... ... @@ -136,13 +136,14 @@ pkginclude_HEADERS = GammaLib.hpp \
136 136 GModelSpatialDiffuseMap.hpp \
137 137 GModelSpectral.hpp \
138 138 GModelSpectralRegistry.hpp \
  139 + GModelSpectralBrokenPlaw.hpp \
  140 + GModelSpectralConst.hpp \
139 141 GModelSpectralExpPlaw.hpp \
140 142 GModelSpectralFunc.hpp \
  143 + GModelSpectralLogParabola.hpp \
141 144 GModelSpectralNodes.hpp \
142 145 GModelSpectralPlaw.hpp \
143 146 GModelSpectralPlaw2.hpp \
144   - GModelSpectralLogParabola.hpp \
145   - GModelSpectralConst.hpp \
146 147 GModelTemporal.hpp \
147 148 GModelTemporalRegistry.hpp \
148 149 GModelTemporalConst.hpp
... ...
inst/lat/test/test_python.py
  View file @4960658
... ... @@ -41,7 +41,7 @@ def analyse_unbinned(xmlname):
41 41 # Load LAT observation
42 42 lat = GLATObservation()
43 43 lat.load_unbinned("data/p7v6/ft1.fits", "data/p7v6/ft2.fits", "data/p7v6/ltcube.fits")
44   - lat.response("P6_v3_diff", "../caldb")
  44 + lat.response("P7REP_SOURCE_V15", "../caldb")
45 45  
46 46 # Setup ROI covered by data
47 47 instDir = GLATInstDir()
... ... @@ -83,7 +83,7 @@ def analyse_binned(xmlname):
83 83 # Load LAT observation
84 84 lat = GLATObservation()
85 85 lat.load_binned("data/p7v6/srcmap.fits", "data/p7v6/binned_expmap.fits", "data/p7v6/ltcube.fits")
86   - lat.response("P6_v3_diff", "../caldb")
  86 + lat.response("P7REP_SOURCE_V15", "../caldb")
87 87 lat.response().force_mean(True)
88 88  
89 89 # Append LAT observation to container
... ... @@ -190,5 +190,5 @@ if __name__ == &#39;__main__&#39;:
190 190 print("**************************************")
191 191  
192 192 # Analyse data
193   - # analyse_unbinned("data/p7v6/crab_model.xml")
194   - analyse_binned("data/p7v6/crab_model.xml")
  193 + analyse_unbinned("data/p7v6/crab_model.xml")
  194 +# analyse_binned("data/p7v6/crab_model.xml")
... ...
pyext/GModelSpectralBrokenPlaw.i 0 → 100644
  View file @4960658
  1 +/***************************************************************************
  2 + * GModelSpectralBrokenPlaw.i - Broken power law spectrum class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Anneli Schulz *
  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 GModelSpectralBrokenPlaw.i
  23 + * @brief Broken power law spectrum class definition
  24 + * @author Anneli Schulz
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "GModelSpectralBrokenPlaw.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class GModelSpectralBrokenPlaw
  34 + *
  35 + * @brief Broken Power law spectral model class
  36 + ***************************************************************************/
  37 +class GModelSpectralBrokenPlaw : public GModelSpectral {
  38 +public:
  39 + // Constructors and destructors
  40 + GModelSpectralBrokenPlaw(void);
  41 + GModelSpectralBrokenPlaw(const double& prefactor,
  42 + const double& index1,
  43 + const GEnergy& breakenergy,
  44 + const double& index2);
  45 + explicit GModelSpectralBrokenPlaw(const GXmlElement& xml);
  46 + GModelSpectralBrokenPlaw(const GModelSpectralBrokenPlaw& model);
  47 + virtual ~GModelSpectralBrokenPlaw(void);
  48 +
  49 + // Implemented pure virtual methods
  50 + virtual void clear(void);
  51 + virtual GModelSpectralBrokenPlaw* clone(void) const;
  52 + virtual std::string type(void) const;
  53 + virtual double eval(const GEnergy& srcEng,
  54 + const GTime& srcTime) const;
  55 + virtual double eval_gradients(const GEnergy& srcEng,
  56 + const GTime& srcTime);
  57 + virtual double flux(const GEnergy& emin,
  58 + const GEnergy& emax) const;
  59 + virtual double eflux(const GEnergy& emin,
  60 + const GEnergy& emax) const;
  61 + virtual GEnergy mc(const GEnergy& emin,
  62 + const GEnergy& emax,
  63 + const GTime& time,
  64 + GRan& ran) const;
  65 + virtual void read(const GXmlElement& xml);
  66 + virtual void write(GXmlElement& xml) const;
  67 +
  68 + // Other methods
  69 + double prefactor(void) const;
  70 + double index1(void) const;
  71 + double index2(void) const;
  72 + GEnergy breakenergy(void) const;
  73 + void prefactor(const double& prefactor);
  74 + void index1(const double& index1);
  75 + void index2(const double& index2);
  76 + void breakenergy(const GEnergy& breakenergy);
  77 +};
  78 +
  79 +
  80 +/***********************************************************************//**
  81 + * @brief GModelSpectralBrokenPlaw class extension
  82 + ***************************************************************************/
  83 +%extend GModelSpectralBrokenPlaw {
  84 + GModelSpectralBrokenPlaw copy() {
  85 + return (*self);
  86 + }
  87 +};
... ...
pyext/gammalib/model.i
  View file @4960658
... ... @@ -73,13 +73,14 @@
73 73 %include "GModelSpatialDiffuseMap.i"
74 74 %include "GModelSpectral.i"
75 75 %include "GModelSpectralRegistry.i"
  76 +%include "GModelSpectralBrokenPlaw.i"
76 77 %include "GModelSpectralConst.i"
77 78 %include "GModelSpectralExpPlaw.i"
78 79 %include "GModelSpectralFunc.i"
  80 +%include "GModelSpectralLogParabola.i"
79 81 %include "GModelSpectralNodes.i"
80 82 %include "GModelSpectralPlaw.i"
81 83 %include "GModelSpectralPlaw2.i"
82   -%include "GModelSpectralLogParabola.i"
83 84 %include "GModelTemporal.i"
84 85 %include "GModelTemporalRegistry.i"
85 86 %include "GModelTemporalConst.i"
... ...
src/model/GModelSpectralBrokenPlaw.cpp 0 → 100644
  View file @4960658
  1 +/***************************************************************************
  2 + * GModelSpectralBrokenPlaw.cpp - Broken power law spectrum class *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2013 by Anneli Schulz *
  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 GModelSpectralBrokenPlaw.cpp
  23 + * @brief Broken power law spectrum class implementation
  24 + * @author Anneli Schulz
  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 "GModelSpectralBrokenPlaw.hpp"
  35 +#include "GModelSpectralRegistry.hpp"
  36 +
  37 +/* __ Constants __________________________________________________________ */
  38 +
  39 +/* __ Globals ____________________________________________________________ */
  40 +const GModelSpectralBrokenPlaw g_spectral_plaw_seed;
  41 +const GModelSpectralRegistry g_spectral_plaw_registry(&g_spectral_plaw_seed);
  42 +
  43 +/* __ Method name definitions ____________________________________________ */
  44 +#define G_MC "GModelSpectralBrokenPlaw::mc(GEnergy&, GEnergy&, GTime&,"\
  45 + " GRan&)"
  46 +#define G_READ "GModelSpectralBrokenPlaw::read(GXmlElement&)"
  47 +#define G_WRITE "GModelSpectralBrokenPlaw::write(GXmlElement&)"
  48 +
  49 +/* __ Macros _____________________________________________________________ */
  50 +
  51 +/* __ Coding definitions _________________________________________________ */
  52 +
  53 +/* __ Debug definitions __________________________________________________ */
  54 +
  55 +
  56 +/*==========================================================================
  57 + = =
  58 + = Constructors/destructors =
  59 + = =
  60 + ==========================================================================*/
  61 +
  62 +/***********************************************************************//**
  63 + * @brief Void constructor
  64 + ***************************************************************************/
  65 +GModelSpectralBrokenPlaw::GModelSpectralBrokenPlaw(void) : GModelSpectral()
  66 +{
  67 + // Initialise private members for clean destruction
  68 + init_members();
  69 +
  70 + // Return
  71 + return;
  72 +}
  73 +
  74 +
  75 +/***********************************************************************//**
  76 + * @brief Parameter constructor
  77 + *
  78 + * @param[in] prefactor Power law pre factor (ph/cm2/s/MeV).
  79 + * @param[in] index1 Power law index1.
  80 + * @param[in] breakenergy break energy.
  81 + * @param[in] index2 Power law index1.
  82 + *
  83 + * Constructs a spectral power law using the model parameters
  84 + *
  85 + * power law @p prefactor (ph/cm2/s/MeV)
  86 + * spectral @p index1
  87 + * Energy @p breakenergy of spectral break
  88 + * spectral @p index2
  89 + ***************************************************************************/
  90 +GModelSpectralBrokenPlaw::GModelSpectralBrokenPlaw(const double& prefactor,
  91 + const double& index1,
  92 + const GEnergy& breakenergy,
  93 + const double& index2) :
  94 + GModelSpectral()
  95 +{
  96 + // Initialise members
  97 + init_members();
  98 +
  99 + // Set parameters
  100 + m_norm.value(prefactor);
  101 + m_index1.value(index1);
  102 + m_breakenergy.value(breakenergy.MeV()); // Internally stored in MeV
  103 + m_index2.value(index2);
  104 +
  105 + // Perform autoscaling of parameter
  106 + autoscale();
  107 +
  108 + // Return
  109 + return;
  110 +}
  111 +
  112 +
  113 +/***********************************************************************//**
  114 + * @brief XML constructor
  115 + *
  116 + * @param[in] xml XML element containing position information.
  117 + *
  118 + * Constructs a power law spectral model by extracting information from an
  119 + * XML element. See the read() method for more information about the expected
  120 + * structure of the XML element.
  121 + ***************************************************************************/
  122 +GModelSpectralBrokenPlaw::GModelSpectralBrokenPlaw(const GXmlElement& xml) :
  123 + GModelSpectral()
  124 +{
  125 + // Initialise members
  126 + init_members();
  127 +
  128 + // Read information from XML element
  129 + read(xml);
  130 +
  131 + // Return
  132 + return;
  133 +}
  134 +
  135 +
  136 +/***********************************************************************//**
  137 + * @brief Copy constructor
  138 + *
  139 + * @param[in] model Broken power law model.
  140 + ***************************************************************************/
  141 +GModelSpectralBrokenPlaw::GModelSpectralBrokenPlaw(const GModelSpectralBrokenPlaw& model) :
  142 + GModelSpectral(model)
  143 +{
  144 + // Initialise members
  145 + init_members();
  146 +
  147 + // Copy members
  148 + copy_members(model);
  149 +
  150 + // Return
  151 + return;
  152 +}
  153 +
  154 +
  155 +/***********************************************************************//**
  156 + * @brief Destructor
  157 + ***************************************************************************/
  158 +GModelSpectralBrokenPlaw::~GModelSpectralBrokenPlaw(void)
  159 +{
  160 + // Free members
  161 + free_members();
  162 +
  163 + // Return
  164 + return;
  165 +}
  166 +
  167 +
  168 +/*==========================================================================
  169 + = =
  170 + = Operators =
  171 + = =
  172 + ==========================================================================*/
  173 +
  174 +/***********************************************************************//**
  175 + * @brief Assignment operator
  176 + *
  177 + * @param[in] model Broken power law model.
  178 + * @return Broken power law model.
  179 + ***************************************************************************/
  180 +GModelSpectralBrokenPlaw& GModelSpectralBrokenPlaw::operator=(const GModelSpectralBrokenPlaw& model)
  181 +{
  182 + // Execute only if object is not identical
  183 + if (this != &model) {
  184 +
  185 + // Copy base class members
  186 + this->GModelSpectral::operator=(model);
  187 +
  188 + // Free members
  189 + free_members();
  190 +
  191 + // Initialise members
  192 + init_members();
  193 +
  194 + // Copy members
  195 + copy_members(model);
  196 +
  197 + } // endif: object was not identical
  198 +
  199 + // Return
  200 + return *this;
  201 +}
  202 +
  203 +
  204 +/*==========================================================================
  205 + = =
  206 + = Public methods =
  207 + = =
  208 + ==========================================================================*/
  209 +
  210 +/***********************************************************************//**
  211 + * @brief Clear broken power law model
  212 + ***************************************************************************/
  213 +void GModelSpectralBrokenPlaw::clear(void)
  214 +{
  215 + // Free class members (base and derived classes, derived class first)
  216 + free_members();
  217 + this->GModelSpectral::free_members();
  218 +
  219 + // Initialise members
  220 + this->GModelSpectral::init_members();
  221 + init_members();
  222 +
  223 + // Return
  224 + return;
  225 +}
  226 +
  227 +
  228 +/***********************************************************************//**
  229 + * @brief Clone broken power law model
  230 + *
  231 + * @return Pointer to deep copy of broken power law model.
  232 + ***************************************************************************/
  233 +GModelSpectralBrokenPlaw* GModelSpectralBrokenPlaw::clone(void) const
  234 +{
  235 + // Clone spectral broken power law model
  236 + return new GModelSpectralBrokenPlaw(*this);
  237 +}
  238 +
  239 +
  240 +/***********************************************************************//**
  241 + * @brief Evaluate function
  242 + *
  243 + * @param[in] srcEng True photon energy.
  244 + * @param[in] srcTime True photon arrival time.
  245 + * @return Model value (ph/cm2/s/MeV).
  246 + *
  247 + * Evaluates
  248 + *
  249 + * \f[
  250 + * S_{\rm E}(E | t) = {\tt m\_norm}
  251 + * \left( \frac{E}{\tt m\_breakenergy} \right)^{\tt m\_index}
  252 + * \f]
  253 + *
  254 + * where
  255 + * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  256 + * - \f${\tt m\_index}\f$ is the spectral index, and
  257 + * - \f${\tt m\_breakenergy}\f$ is the breakenergy energy.
  258 + *
  259 + * @todo The method expects that energy!=0. Otherwise Inf or NaN may result.
  260 + ***************************************************************************/
  261 +double GModelSpectralBrokenPlaw::eval(const GEnergy& srcEng,
  262 + const GTime& srcTime) const
  263 +{
  264 + // Update the evaluation cache
  265 + update_eval_cache(srcEng);
  266 +
  267 + // Compute function value
  268 + double value = m_norm.value() * m_last_power;
  269 +
  270 + // Compile option: Check for NaN/Inf
  271 + #if defined(G_NAN_CHECK)
  272 + if (gammalib::isnotanumber(value) || gammalib::isinfinite(value)) {
  273 + std::cout << "*** ERROR: GModelSpectralBrokenPlaw::eval";
  274 + std::cout << "(srcEng=" << srcEng;
  275 + std::cout << ", srcTime=" << srcTime << "):";
  276 + std::cout << " NaN/Inf encountered";
  277 + std::cout << " (value=" << value;
  278 + std::cout << ", m_norm=" << m_norm.value();
  279 + std::cout << ", m_index1=" << m_index1.value();
  280 + std::cout << ", m_breakenergy=" << m_breakenergy.value();
  281 + std::cout << ", m_index2=" << m_index2.value();
  282 + std::cout << ", m_last_power=" << m_last_power;
  283 + std::cout << ")" << std::endl;
  284 + }
  285 + #endif
  286 +
  287 + // Return
  288 + return value;
  289 +}
  290 +
  291 +
  292 +/***********************************************************************//**
  293 + * @brief Evaluate function and gradients
  294 + *
  295 + * @param[in] srcEng True photon energy.
  296 + * @param[in] srcTime True photon arrival time.
  297 + * @return Model value (ph/cm2/s/MeV).
  298 + *
  299 + * Evaluates
  300 + *
  301 + * \f[
  302 + * S_{\rm E}(E | t) = {\tt m\_norm}
  303 + * \left( \frac{E}{\tt m\_breakenergy} \right)^{\tt m\_index}
  304 + * \f]
  305 + *
  306 + * where
  307 + * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
  308 + * - \f${\tt m\_index}\f$ is the spectral index, and
  309 + * - \f${\tt m\_breakenergy}\f$ is the breakenergy energy.
  310 + *
  311 + * The method also evaluates the partial derivatives of the model with
  312 + * respect to the parameters using
  313 + *
  314 + * \f[
  315 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} =
  316 + * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}}
  317 + * \f]
  318 + *
  319 + * \f[
  320 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} =
  321 + * S_{\rm E}(E | t) \, \ln(E/{\tt m_breakenergy})
  322 + * \f]
  323 + *
  324 + * \f[
  325 + * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_breakenergy}} =
  326 + * -S_{\rm E}(E | t) \,
  327 + * \left( \frac{{\tt m\_index}}{{\tt m\_breakenergy}} \right)
  328 + * \f]
  329 + *
  330 + * @todo The method expects that energy!=0. Otherwise Inf or NaN may result.
  331 + ***************************************************************************/
  332 +double GModelSpectralBrokenPlaw::eval_gradients(const GEnergy& srcEng,
  333 + const GTime& srcTime)
  334 +{
  335 + // Update the evaluation cache
  336 + update_eval_cache(srcEng);
  337 +
  338 + // Compute function value
  339 + double value = m_norm.value() * m_last_power;
  340 +
  341 + // Compute normalisation gradient
  342 + double g_norm = (m_norm.isfree())
  343 + ? m_norm.scale() * m_last_power : 0.0;
  344 + m_norm.factor_gradient(g_norm);
  345 +
  346 + // Compute index and break value gradients
  347 + if (srcEng.MeV() < m_breakenergy.value()) {
  348 + double g_index = (m_index1.isfree())
  349 + ? value * m_index1.scale() * m_last_log_e_norm
  350 + : 0.0;
  351 + double g_break = (m_breakenergy.isfree())
  352 + ? -value * m_last_index1 / m_breakenergy.factor_value()
  353 + : 0.0;
  354 + m_index1.factor_gradient(g_index);
  355 + m_index2.factor_gradient(0.0);
  356 + m_breakenergy.factor_gradient(g_break);
  357 + }
  358 + else {
  359 + double g_index = (m_index2.isfree())
  360 + ? value * m_index2.scale() * m_last_log_e_norm
  361 + : 0.0;
  362 + double g_break = (m_breakenergy.isfree())
  363 + ? -value * m_last_index2 / m_breakenergy.factor_value()
  364 + : 0.0;
  365 + m_index1.factor_gradient(0.0);
  366 + m_index2.factor_gradient(g_index);
  367 + m_breakenergy.factor_gradient(g_break);
  368 + }
  369 +
  370 + // Compile option: Check for NaN/Inf
  371 + #if defined(G_NAN_CHECK)
  372 + if (gammalib::isnotanumber(value) || gammalib::isinfinite(value)) {
  373 + std::cout << "*** ERROR: GModelSpectralBrokenPlaw::eval_gradients";
  374 + std::cout << "(srcEng=" << srcEng;
  375 + std::cout << ", srcTime=" << srcTime << "):";
  376 + std::cout << " NaN/Inf encountered";
  377 + std::cout << " (value=" << value;
  378 + std::cout << ", e_norm=" << m_last_e_norm;
  379 + std::cout << ", power=" << m_last_power;
  380 + std::cout << ", g_norm=" << g_norm;
  381 + std::cout << ")" << std::endl;
  382 + }
  383 + #endif
  384 +
  385 + // Return
  386 + return value;
  387 +}
  388 +
  389 +
  390 +/***********************************************************************//**
  391 + * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
  392 + *
  393 + * @param[in] emin Minimum photon energy.
  394 + * @param[in] emax Maximum photon energy.
  395 + * @return Photon flux (ph/cm2/s).
  396 + *
  397 + * Computes
  398 + *
  399 + * \f[
  400 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
  401 + * \f]
  402 + *
  403 + * where
  404 + * - [@p emin, @p emax] is an energy interval, and
  405 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  406 + * The integration is done analytically.
  407 + ***************************************************************************/
  408 +double GModelSpectralBrokenPlaw::flux(const GEnergy& emin,
  409 + const GEnergy& emax) const
  410 +{
  411 + // Initialise flux
  412 + double flux = 0.0;
  413 +
  414 + // Compute only if integration range is valid
  415 + if (emin < emax) {
  416 +
  417 + // First case: emin < breakenergy < emax
  418 + if (emin.MeV() < m_breakenergy.value() && m_breakenergy.value() < emax.MeV()) {
  419 +
  420 + // Compute photon flux
  421 + flux = m_norm.value() *
  422 + (gammalib::plaw_photon_flux(emin.MeV(),
  423 + m_breakenergy.value(),
  424 + m_breakenergy.value(),
  425 + m_index1.value()) +
  426 + gammalib::plaw_photon_flux(m_breakenergy.value(),
  427 + emax.MeV(),
  428 + m_breakenergy.value(),
  429 + m_index2.value()));
  430 + }
  431 +
  432 + // Second case: breakenergy > emax:
  433 + else if (m_breakenergy.value() > emax.MeV()) {
  434 + flux = m_norm.value() *
  435 + gammalib::plaw_photon_flux(emin.MeV(),
  436 + emax.MeV(),
  437 + m_breakenergy.value(),
  438 + m_index1.value());
  439 + }
  440 +
  441 + // Third case breakenergy < emin:
  442 + else if (m_breakenergy.value() < emin.MeV()) {
  443 + flux = m_norm.value() *
  444 + gammalib::plaw_photon_flux(emin.MeV(),
  445 + emax.MeV(),
  446 + m_breakenergy.value(),
  447 + m_index2.value());
  448 + }
  449 +
  450 + } // endif: integration range was valid
  451 +
  452 + // Return flux
  453 + return flux;
  454 +}
  455 +
  456 +
  457 +/***********************************************************************//**
  458 + * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
  459 + *
  460 + * @param[in] emin Minimum photon energy.
  461 + * @param[in] emax Maximum photon energy.
  462 + * @return Energy flux (erg/cm2/s).
  463 + *
  464 + * Computes
  465 + *
  466 + * \f[
  467 + * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
  468 + * \f]
  469 + *
  470 + * where
  471 + * - [@p emin, @p emax] is an energy interval, and
  472 + * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
  473 + * The integration is done analytically.
  474 + ***************************************************************************/
  475 +double GModelSpectralBrokenPlaw::eflux(const GEnergy& emin,
  476 + const GEnergy& emax) const
  477 +{
  478 + // Initialise flux
  479 + double eflux = 0.0;
  480 +
  481 + // Compute only if integration range is valid
  482 + if (emin < emax) {
  483 +
  484 + // First case: emin < breakenergy < emax
  485 + if (emin.MeV() < m_breakenergy.value() && m_breakenergy.value() < emax.MeV()) {
  486 + eflux = m_norm.value() *
  487 + (gammalib::plaw_energy_flux(emin.MeV(),
  488 + m_breakenergy.value(),
  489 + m_breakenergy.value(),
  490 + m_index1.value()) +
  491 + gammalib::plaw_energy_flux(m_breakenergy.value(),
  492 + emax.MeV(),
  493 + m_breakenergy.value(),
  494 + m_index2.value()));
  495 + }
  496 +
  497 + // Second case: breakenergy > emax:
  498 + else if (m_breakenergy.value() > emax.MeV()) {
  499 + eflux = m_norm.value() *
  500 + gammalib::plaw_energy_flux(emin.MeV(),
  501 + emax.MeV(),
  502 + m_breakenergy.value(),
  503 + m_index1.value());
  504 + }
  505 +
  506 + // Third case breakenergy < emin:
  507 + else if (m_breakenergy.value() < emin.MeV()) {
  508 + eflux = m_norm.value() *
  509 + gammalib::plaw_energy_flux(emin.MeV(),
  510 + emax.MeV(),
  511 + m_breakenergy.value(),
  512 + m_index2.value());
  513 + }
  514 + eflux *= gammalib::MeV2erg;
  515 +
  516 + } // endif: integration range was valid
  517 +
  518 + // Return flux
  519 + return eflux;
  520 +}
  521 +
  522 +
  523 +/***********************************************************************//**
  524 + * @brief Returns Monte Carlo energy between [emin, emax]
  525 + *
  526 + * @param[in] emin Minimum photon energy.
  527 + * @param[in] emax Maximum photon energy.
  528 + * @param[in] time True photon arrival time.
  529 + * @param[in,out] ran Random number generator.
  530 + * @return Energy.
  531 + *
  532 + * @exception GException::erange_invalid
  533 + * Energy range is invalid (emin < emax required).
  534 + *
  535 + * Returns Monte Carlo energy by randomly drawing from a broken power law.
  536 + ***************************************************************************/
  537 +GEnergy GModelSpectralBrokenPlaw::mc(const GEnergy& emin,
  538 + const GEnergy& emax,
  539 + const GTime& time,
  540 + GRan& ran) const
  541 +{
  542 + // Throw an exception if energy range is invalid
  543 + if (emin >= emax) {
  544 + throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
  545 + "Minimum energy < maximum energy required.");
  546 + }
  547 +
  548 + // Allocate energy
  549 + GEnergy energy;
  550 +
  551 + // Update cache
  552 + update_mc_cache(emin, emax);
  553 +
  554 + // Determine in which bin we reside
  555 + int inx = 0;
  556 + if (m_mc_cum.size() > 1) {
  557 + double u = ran.uniform();
  558 + for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
  559 + if (m_mc_cum[inx-1] <= u)
  560 + break;
  561 + }
  562 + }
  563 +
  564 + // Get random energy for specific bin
  565 + if (m_mc_exp[inx] != 0.0) {
  566 + double e_min = m_mc_min[inx];
  567 + double e_max = m_mc_max[inx];
  568 + double u = ran.uniform();
  569 + double eng = (u > 0.0)
  570 + ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
  571 + : 0.0;
  572 + energy.MeV(eng);
  573 + }
  574 + else {
  575 + double e_min = m_mc_min[inx];
  576 + double e_max = m_mc_max[inx];
  577 + double u = ran.uniform();
  578 + double eng = std::exp(u * (e_max - e_min) + e_min);
  579 + energy.MeV(eng);
  580 + }
  581 +
  582 + // Return energy
  583 + return energy;
  584 +}
  585 +
  586 +
  587 +/***********************************************************************//**
  588 + * @brief Read model from XML element
  589 + *
  590 + * @param[in] xml XML element containing power law model information.
  591 + *
  592 + * @exception GException::model_invalid_parnum
  593 + * Invalid number of model parameters found in XML element.
  594 + * @exception GException::model_invalid_parnames
  595 + * Invalid model parameter names found in XML element.
  596 + *
  597 + * Reads the spectral information from an XML element. The format of the XML
  598 + * elements is
  599 + *
  600 + * <spectrum type="BrokenPowerLaw">
  601 + * <parameter name="Prefactor" scale=".." value=".." max=".." min=".." free=".."/>
  602 + * <parameter name="Index1" scale=".." value=".." max=".." min=".." free=".."/>
  603 + * <parameter name="BreakValue" scale=".." value=".." max=".." min=".." free=".."/>
  604 + * <parameter name="Index2" scale=".." value=".." max=".." min=".." free=".."/>
  605 + * </spectrum>
  606 + *
  607 + * @todo Add parameter validity check
  608 + ***************************************************************************/
  609 +void GModelSpectralBrokenPlaw::read(const GXmlElement& xml)
  610 +{
  611 + // Verify that XML element has exactly 4 parameters
  612 + if (xml.elements() != 4 || xml.elements("parameter") != 4) {
  613 + throw GException::model_invalid_parnum(G_READ, xml,
  614 + "Broken Power law model requires exactly 4 parameters.");
  615 + }
  616 +
  617 + // Extract model parameters
  618 + int npar[] = {0, 0, 0, 0};
  619 + for (int i = 0; i < 4; ++i) {
  620 +
  621 + // Get parameter element
  622 + const GXmlElement* par = xml.element("parameter", i);
  623 +
  624 + // Handle prefactor
  625 + if (par->attribute("name") == "Prefactor") {
  626 + m_norm.read(*par);
  627 + npar[0]++;
  628 + }
  629 +
  630 + // Handle index1
  631 + else if (par->attribute("name") == "Index1") {
  632 + m_index1.read(*par);
  633 + npar[1]++;
  634 + }
  635 +
  636 + // Handle breakenergy energy
  637 + else if (par->attribute("name") == "BreakValue") {
  638 + m_breakenergy.read(*par);
  639 + npar[2]++;
  640 + }
  641 + // Handle index2
  642 + else if (par->attribute("name") == "Index2") {
  643 + m_index2.read(*par);
  644 + npar[3]++;
  645 + }
  646 +
  647 +
  648 + } // endfor: looped over all parameters
  649 +
  650 + // Verify that all parameters were found
  651 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] !=1) {
  652 + throw GException::model_invalid_parnames(G_READ, xml,
  653 + "Require \"Prefactor\", \"Index1\", \"BreakValue\" and \"Index2\" parameters.");
  654 + }
  655 +
  656 + // Return
  657 + return;
  658 +}
  659 +
  660 +
  661 +/***********************************************************************//**
  662 + * @brief Write model into XML element
  663 + *
  664 + * @param[in] xml XML element into which model information is written.
  665 + *
  666 + * @exception GException::model_invalid_
  667 + * Existing XML element is not of type "BrokenPowerLaw"
  668 + * @exception GException::model_invalid_parnum
  669 + * Invalid number of model parameters or nodes found in XML element.
  670 + * @exception GException::model_invalid_parnames
  671 + * Invalid model parameter names found in XML element.
  672 + *
  673 + * Writes the spectral information into an XML element. The format of the XML
  674 + * element is
  675 + * * <spectrum type="BrokenPowerLaw">
  676 + * <parameter name="Prefactor" scale=".." value=".." max=".." min=".." free=".."/>
  677 + * <parameter name="Index1" scale=".." value=".." max=".." min=".." free=".."/>
  678 + * <parameter name="BreakValue" scale=".." value=".." max=".." min=".." free=".."/>
  679 + * <parameter name="Index2" scale=".." value=".." max=".." min=".." free=".."/>
  680 + * </spectrum>
  681 + ***************************************************************************/
  682 +void GModelSpectralBrokenPlaw::write(GXmlElement& xml) const
  683 +{
  684 + // Set model type
  685 + if (xml.attribute("type") == "") {
  686 + xml.attribute("type", "BrokenPowerLaw");
  687 + }
  688 +
  689 + // Verify model type
  690 + if (xml.attribute("type") != "BrokenPowerLaw") {
  691 + throw GException::model_invalid_spectral(G_WRITE, xml.attribute("type"),
  692 + "Spectral model is not of type \"BrokenPowerLaw\".");
  693 + }
  694 +
  695 + // If XML element has 0 nodes then append 4 parameter nodes
  696 + if (xml.elements() == 0) {
  697 + xml.append(GXmlElement("parameter name=\"Prefactor\""));
  698 + xml.append(GXmlElement("parameter name=\"Index1\""));
  699 + xml.append(GXmlElement("parameter name=\"BreakValue\""));
  700 + xml.append(GXmlElement("parameter name=\"Index2\""));
  701 + }
  702 +
  703 + // Verify that XML element has exactly 4 parameters
  704 + if (xml.elements() != 4 || xml.elements("parameter") != 4) {
  705 + throw GException::model_invalid_parnum(G_WRITE, xml,
  706 + "Broken Power law model requires exactly 4 parameters.");
  707 + }
  708 +
  709 + // Set or update model parameter attributes
  710 + int npar[] = {0, 0, 0, 0};
  711 + for (int i = 0; i < 4; ++i) {
  712 +
  713 + // Get parameter element
  714 + GXmlElement* par = xml.element("parameter", i);
  715 +
  716 + // Handle prefactor
  717 + if (par->attribute("name") == "Prefactor") {
  718 + npar[0]++;
  719 + m_norm.write(*par);
  720 + }
  721 +
  722 + // Handle index1
  723 + else if (par->attribute("name") == "Index1") {
  724 + npar[1]++;
  725 + m_index1.write(*par);
  726 + }
  727 +
  728 + // Handle break energy
  729 + else if (par->attribute("name") == "BreakValue") {
  730 + npar[2]++;
  731 + m_breakenergy.write(*par);
  732 + }
  733 +
  734 + // Handle index2
  735 + else if (par->attribute("name") == "Index2") {
  736 + npar[3]++;
  737 + m_index2.write(*par);
  738 + }
  739 +
  740 + } // endfor: looped over all parameters
  741 +
  742 + // Check of all required parameters are present
  743 + if (npar[0] != 1 || npar[1] != 1 || npar[2] != 1 || npar[3] != 1) {
  744 + throw GException::model_invalid_parnames(G_WRITE, xml,
  745 + "Require \"Prefactor\", \"Index1\", \"BreakValue\" and \"Index2\" parameters.");
  746 + }
  747 +
  748 + // Return
  749 + return;
  750 +}
  751 +
  752 +
  753 +/***********************************************************************//**
  754 + * @brief Print brokenpowerlaw information
  755 + *
  756 + * @param[in] chatter Chattiness (defaults to NORMAL).
  757 + * @return String containing model information.
  758 + ***************************************************************************/
  759 +std::string GModelSpectralBrokenPlaw::print(const GChatter& chatter) const
  760 +{
  761 + // Initialise result string
  762 + std::string result;
  763 +
  764 + // Continue only if chatter is not silent
  765 + if (chatter != SILENT) {
  766 +
  767 + // Append header
  768 + result.append("=== GModelSpectralBrokenPlaw ===");
  769 +
  770 + // Append information
  771 + result.append("\n"+gammalib::parformat("Number of parameters"));
  772 + result.append(gammalib::str(size()));
  773 + for (int i = 0; i < size(); ++i) {
  774 + result.append("\n"+m_pars[i]->print(chatter));
  775 + }
  776 +
  777 + } // endif: chatter was not silent
  778 +
  779 + // Return result
  780 + return result;
  781 +}
  782 +
  783 +
  784 +/*==========================================================================
  785 + = =
  786 + = Private methods =
  787 + = =
  788 + ==========================================================================*/
  789 +
  790 +/***********************************************************************//**
  791 + * @brief Initialise class members
  792 + ***************************************************************************/
  793 +void GModelSpectralBrokenPlaw::init_members(void)
  794 +{
  795 + // Initialise powerlaw normalisation
  796 + m_norm.clear();
  797 + m_norm.name("Prefactor");
  798 + m_norm.unit("ph/cm2/s/MeV");
  799 + m_norm.scale(1.0);
  800 + m_norm.value(1.0); // default: 1.0
  801 + m_norm.min(0.0); // min: 0.0
  802 + m_norm.free();
  803 + m_norm.gradient(0.0);
  804 + m_norm.hasgrad(true);
  805 +
  806 + // Initialise powerlaw index1
  807 + m_index1.clear();
  808 + m_index1.name("Index1");
  809 + m_index1.scale(1.0);
  810 + m_index1.value(-2.0); // default: -2.0
  811 + m_index1.range(-10.0,+10.0); // range: [-10,+10]
  812 + m_index1.free();
  813 + m_index1.gradient(0.0);
  814 + m_index1.hasgrad(true);
  815 +
  816 + // Initialise powerlaw index2
  817 + m_index2.clear();
  818 + m_index2.name("Index2");
  819 + m_index2.scale(1.0);
  820 + m_index2.value(-2.0); // default: -2.0
  821 + m_index2.range(-10.0,+10.0); // range: [-10,+10]
  822 + m_index2.free();
  823 + m_index2.gradient(0.0);
  824 + m_index2.hasgrad(true);
  825 +
  826 + // Initialise break energy
  827 + m_breakenergy.clear();
  828 + m_breakenergy.name("BreakValue");
  829 + m_breakenergy.unit("MeV");
  830 + m_breakenergy.scale(1.0);
  831 + m_breakenergy.value(100.0); // default: 100
  832 + m_breakenergy.fix();
  833 + m_breakenergy.gradient(0.0);
  834 + m_breakenergy.hasgrad(true);
  835 +
  836 + // Set parameter pointer(s)
  837 + m_pars.clear();
  838 + m_pars.push_back(&m_norm);
  839 + m_pars.push_back(&m_index1);
  840 + m_pars.push_back(&m_breakenergy);
  841 + m_pars.push_back(&m_index2);
  842 +
  843 + // Initialise eval cache
  844 + m_last_energy.clear();
  845 + m_last_index1 = 1.0e30;
  846 + m_last_index2 = 1.0e30;
  847 + m_last_breakenergy = 1.0e30;
  848 + m_last_e_norm = 0.0;
  849 + m_last_log_e_norm = 0.0;
  850 + m_last_power = 0.0;
  851 +
  852 + // Initialise MC cache
  853 + m_mc_exponent1 = 0.0;
  854 + m_mc_exponent2 = 0.0;
  855 + m_mc_pow_emin = 0.0;
  856 + m_mc_pow_ewidth = 0.0;
  857 +
  858 + // Initialise MC cache
  859 + m_mc_emin = 0.0;
  860 + m_mc_emax = 0.0;
  861 +
  862 +
  863 + // Return
  864 + return;
  865 +}
  866 +
  867 +
  868 +/***********************************************************************//**
  869 + * @brief Copy class members
  870 + *
  871 + * @param[in] model GModelSpectralBrokenPlaw members which should be copied.
  872 + ***************************************************************************/
  873 +void GModelSpectralBrokenPlaw::copy_members(const GModelSpectralBrokenPlaw& model)
  874 +{
  875 + // Copy members
  876 + m_norm = model.m_norm;
  877 + m_index1 = model.m_index1;
  878 + m_breakenergy = model.m_breakenergy;
  879 + m_index2 = model.m_index2;
  880 +
  881 + // Set parameter pointer(s)
  882 + m_pars.clear();
  883 + m_pars.push_back(&m_norm);
  884 + m_pars.push_back(&m_index1);
  885 + m_pars.push_back(&m_breakenergy);
  886 + m_pars.push_back(&m_index2);
  887 +
  888 + // Copy eval cache
  889 + m_last_energy = model.m_last_energy;
  890 + m_last_index1 = model.m_last_index1;
  891 + m_last_breakenergy = model.m_last_breakenergy;
  892 + m_last_index2 = model.m_last_index2;
  893 + m_last_e_norm = model.m_last_e_norm;
  894 + m_last_log_e_norm = model.m_last_log_e_norm;
  895 + m_last_power = model.m_last_power;
  896 +
  897 + // Copy MC cache
  898 + m_mc_emin = model.m_mc_emin;
  899 + m_mc_emax = model.m_mc_emax;
  900 + m_mc_exponent1 = model.m_mc_exponent1;
  901 + m_mc_exponent2 = model.m_mc_exponent2;
  902 + m_mc_pow_emin = model.m_mc_pow_emin;
  903 + m_mc_pow_ewidth = model.m_mc_pow_ewidth;
  904 +
  905 + // Return
  906 + return;
  907 +}
  908 +
  909 +
  910 +/***********************************************************************//**
  911 + * @brief Delete class members
  912 + ***************************************************************************/
  913 +void GModelSpectralBrokenPlaw::free_members(void)
  914 +{
  915 + // Return
  916 + return;
  917 +}
  918 +
  919 +
  920 +/***********************************************************************//**
  921 + * @brief Update eval precomputation cache
  922 + *
  923 + * @param[in] energy Energy.
  924 + *
  925 + * Updates the precomputation cache for eval() and eval_gradients() methods.
  926 + ***************************************************************************/
  927 +void GModelSpectralBrokenPlaw::update_eval_cache(const GEnergy& energy) const
  928 +{
  929 + // Get parameter values (takes 2 multiplications which are difficult
  930 + // to avoid)
  931 + double index1 = m_index1.value();
  932 + double breakenergy = m_breakenergy.value();
  933 + double index2 = m_index2.value();
  934 +
  935 + // If the energy or one of the parameters index1, index2 or breakenergy
  936 + // energy has changed then recompute the cache
  937 + if ((m_last_energy != energy) ||
  938 + (m_last_index1 != index1) ||
  939 + (m_last_index2 != index2) ||
  940 + (m_last_breakenergy != breakenergy)) {
  941 +
  942 + // Store actual energy and parameter values
  943 + m_last_energy = energy;
  944 + m_last_index1 = index1;
  945 + m_last_index2 = index2;
  946 + m_last_breakenergy = breakenergy;
  947 +
  948 + // Compute and store value
  949 + double eng = energy.MeV();
  950 + m_last_e_norm = eng / m_last_breakenergy;
  951 + m_last_log_e_norm = std::log(m_last_e_norm);
  952 + if (eng < m_last_breakenergy ) {
  953 + m_last_power = std::pow(m_last_e_norm, m_last_index1);
  954 + }
  955 + else {
  956 + m_last_power = std::pow(m_last_e_norm, m_last_index2);
  957 + }
  958 +
  959 + } // endif: recomputation was required
  960 +
  961 + // Return
  962 + return;
  963 +}
  964 +
  965 +
  966 +/***********************************************************************//**
  967 + * @brief Update Monte Carlo pre computation cache
  968 + *
  969 + * @param[in] emin Minimum photon energy.
  970 + * @param[in] emax Maximum photon energy.
  971 + *
  972 + * Updates the precomputation cache for Monte Carlo simulations.
  973 + ***************************************************************************/
  974 +void GModelSpectralBrokenPlaw::update_mc_cache(const GEnergy& emin,
  975 + const GEnergy& emax) const
  976 +{
  977 + // Check if we need to update the cache
  978 + if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) {
  979 +
  980 + // Store new energy interval
  981 + m_mc_emin = emin.MeV();
  982 + m_mc_emax = emax.MeV();
  983 +
  984 + // Initialise cache
  985 + m_mc_cum.clear();
  986 + m_mc_min.clear();
  987 + m_mc_max.clear();
  988 + m_mc_exp.clear();
  989 +
  990 + // Get energy range in MeV
  991 + double e_min = emin.MeV();
  992 + double e_max = emax.MeV();
  993 +
  994 + // Continue only if e_max > e_min
  995 + if (e_max > e_min) {
  996 +
  997 + // Allocate flux
  998 + double flux;
  999 +
  1000 + // Determine left node index for minimum and maximum energy
  1001 + int inx_emin = (e_min < m_breakenergy.value() ) ? 0 : 1;
  1002 + int inx_emax = (e_max < m_breakenergy.value() ) ? 0 : 1;
  1003 +
  1004 + // If both energies are within the same node then just
  1005 + // add this one node on the stack
  1006 + if (inx_emin == inx_emax) {
  1007 + double exp_valid = (e_min < m_breakenergy.value())
  1008 + ? m_index1.value() : m_index2.value();
  1009 + flux = m_norm.value() *
  1010 + gammalib::plaw_photon_flux(e_min,
  1011 + e_max,
  1012 + m_breakenergy.value(),
  1013 + exp_valid);
  1014 + m_mc_cum.push_back(flux);
  1015 + m_mc_min.push_back(e_min);
  1016 + m_mc_max.push_back(e_max);
  1017 + m_mc_exp.push_back(exp_valid);
  1018 + }
  1019 +
  1020 + // ... otherwise integrate over both nodes
  1021 + else {
  1022 + // just enter the values for first pl: bin [0]
  1023 + flux = m_norm.value() *
  1024 + gammalib::plaw_photon_flux(e_min,
  1025 + m_breakenergy.value(),
  1026 + m_breakenergy.value(),
  1027 + m_index1.value());
  1028 + m_mc_cum.push_back(flux);
  1029 + m_mc_exp.push_back(m_index1.value());
  1030 + m_mc_min.push_back(e_min);
  1031 + m_mc_max.push_back(m_breakenergy.value());
  1032 +
  1033 + // and for
  1034 + flux = m_norm.value() *
  1035 + gammalib::plaw_photon_flux(m_breakenergy.value(),
  1036 + e_max,
  1037 + m_breakenergy.value(),
  1038 + m_index2.value());
  1039 + m_mc_cum.push_back(flux);
  1040 + m_mc_exp.push_back(m_index2.value());
  1041 + m_mc_max.push_back(e_max);
  1042 + m_mc_min.push_back(m_breakenergy.value());
  1043 + } // endelse: emin and emax not between same nodes
  1044 +
  1045 + // Build cumulative distribution
  1046 + for (int i = 1; i < m_mc_cum.size(); ++i) {
  1047 + m_mc_cum[i] += m_mc_cum[i-1];
  1048 + }
  1049 + double norm = m_mc_cum[m_mc_cum.size()-1];
  1050 + for (int i = 0; i < m_mc_cum.size(); ++i) {
  1051 + m_mc_cum[i] /= norm;
  1052 + }
  1053 +
  1054 + // Set MC values
  1055 + for (int i = 0; i < m_mc_cum.size(); ++i) {
  1056 +
  1057 + // Compute exponent
  1058 + double exponent = m_mc_exp[i] + 1.0;
  1059 +
  1060 + // Exponent dependend computation
  1061 + if (std::abs(exponent) > 1.0e-11) {
  1062 + m_mc_exp[i] = exponent;
  1063 + m_mc_min[i] = std::pow(m_mc_min[i], exponent);
  1064 + m_mc_max[i] = std::pow(m_mc_max[i], exponent);
  1065 + }
  1066 + else {
  1067 + m_mc_exp[i] = 0.0;
  1068 + m_mc_min[i] = std::log(m_mc_min[i]);
  1069 + m_mc_max[i] = std::log(m_mc_max[i]);
  1070 + }
  1071 +
  1072 + } // endfor: set MC values
  1073 +
  1074 + } // endif: e_max > e_min
  1075 +
  1076 + } // endif: Update was required
  1077 +
  1078 + // Return
  1079 + return;
  1080 +}
... ...
src/model/Makefile.am
  View file @4960658
... ... @@ -27,13 +27,14 @@ sources = GModelPar.cpp \
27 27 GModelSpatialDiffuseMap.cpp \
28 28 GModelSpectral.cpp \
29 29 GModelSpectralRegistry.cpp \
  30 + GModelSpectralBrokenPlaw.cpp \
30 31 GModelSpectralConst.cpp \
31 32 GModelSpectralExpPlaw.cpp \
32 33 GModelSpectralFunc.cpp \
  34 + GModelSpectralLogParabola.cpp \
33 35 GModelSpectralNodes.cpp \
34 36 GModelSpectralPlaw.cpp \
35 37 GModelSpectralPlaw2.cpp \
36   - GModelSpectralLogParabola.cpp \
37 38 GModelTemporal.cpp \
38 39 GModelTemporalRegistry.cpp \
39 40 GModelTemporalConst.cpp \
... ...
test/data/model_point_bplaw.xml 0 → 100644
  View file @4960658
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Crab" type="PointSource">
  4 + <spectrum type="BrokenPowerLaw">
  5 + <parameter name="Prefactor" scale="1e-16" value="5.7" min="1e-07" max="1000.0" free="1"/>
  6 + <parameter name="Index1" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
  7 + <parameter name="BreakValue" scale="1e6" value="0.3" min="0.01" max="1000.0" free="1"/>
  8 + <parameter name="Index2" scale="-1" value="2.70" min="0.01" max="+5.0" free="1"/>
  9 + </spectrum>
  10 + <spatialModel type="SkyDirFunction">
  11 + <parameter free="0" max="360" min="-360" name="RA" scale="1" value="83.6331" />
  12 + <parameter free="0" max="90" min="-90" name="DEC" scale="1" value="22.0145" />
  13 + </spatialModel>
  14 + </source>
  15 +</source_library>
... ...
test/test_GModel.cpp
  View file @4960658
... ... @@ -52,6 +52,7 @@ void TestGModel::set(void)
52 52 m_xml_model_point_plaw = "data/model_point_plaw.xml";
53 53 m_xml_model_point_plaw2 = "data/model_point_plaw2.xml";
54 54 m_xml_model_point_eplaw = "data/model_point_eplaw.xml";
  55 + m_xml_model_point_bplaw = "data/model_point_bplaw.xml";
55 56 m_xml_model_point_logparabola = "data/model_point_logparabola.xml";
56 57 m_xml_model_point_nodes = "data/model_point_nodes.xml";
57 58 m_xml_model_point_filefct = "data/model_point_filefct.xml";
... ... @@ -83,6 +84,7 @@ void TestGModel::set(void)
83 84 append(static_cast<pfunction>(&TestGModel::test_plaw), "Test GModelSpectralPlaw");
84 85 append(static_cast<pfunction>(&TestGModel::test_plaw2), "Test GModelSpectralPlaw2");
85 86 append(static_cast<pfunction>(&TestGModel::test_eplaw), "Test GModelSpectralExpPlaw");
  87 + append(static_cast<pfunction>(&TestGModel::test_bplaw), "Test GModelSpectralBrokenPlaw");
86 88 append(static_cast<pfunction>(&TestGModel::test_logparabola), "Test GModelSpectralLogParabola");
87 89 append(static_cast<pfunction>(&TestGModel::test_nodes), "Test GModelSpectralNodes");
88 90 append(static_cast<pfunction>(&TestGModel::test_filefct), "Test GModelSpectralFunc");
... ... @@ -1406,6 +1408,92 @@ void TestGModel::test_eplaw(void)
1406 1408  
1407 1409  
1408 1410 /***********************************************************************//**
  1411 + * @brief Test GModelSpectralBrokenPlaw class
  1412 + ***************************************************************************/
  1413 +void TestGModel::test_bplaw(void)
  1414 +{
  1415 + // Test void constructor
  1416 + test_try("Test void constructor");
  1417 + try {
  1418 + GModelSpectralBrokenPlaw model;
  1419 + test_assert(model.type() == "BrokenPowerLaw",
  1420 + "Model type \"BrokenPowerLaw\" expected.");
  1421 + test_try_success();
  1422 + }
  1423 + catch (std::exception &e) {
  1424 + test_try_failure(e);
  1425 + }
  1426 +
  1427 + // Test value constructor
  1428 + test_try("Test value constructor");
  1429 + try {
  1430 + GModelSpectralBrokenPlaw model(2.0, -2.1, GEnergy(100.0, "MeV"), -2.8);
  1431 + test_value(model.prefactor(), 2.0);
  1432 + test_value(model.index1(), -2.1);
  1433 + test_value(model.breakenergy().MeV(), 100.0);
  1434 + test_value(model.index2(), -2.8);
  1435 + test_try_success();
  1436 + }
  1437 + catch (std::exception &e) {
  1438 + test_try_failure(e);
  1439 + }
  1440 +
  1441 + // Test XML constructor and value
  1442 + test_try("Test XML constructor, value and gradients");
  1443 + try {
  1444 + // Test XML constructor
  1445 + GXml xml(m_xml_model_point_bplaw);
  1446 + GXmlElement* element = xml.element(0)->element(0)->element("spectrum", 0);
  1447 + GModelSpectralBrokenPlaw model(*element);
  1448 + test_value(model.size(), 4);
  1449 + test_assert(model.type() == "BrokenPowerLaw", "Expected \"BrokenPowerLaw\"");
  1450 + test_value(model.prefactor(), 5.7e-16);
  1451 + test_value(model.index1(), -2.48);
  1452 + test_value(model.breakenergy().TeV(), 0.3);
  1453 + test_value(model.index2(), -2.70);
  1454 +
  1455 + // Test prefactor method
  1456 + model.prefactor(2.3e-16);
  1457 + test_value(model.prefactor(), 2.3e-16);
  1458 +
  1459 + // Test index1 method
  1460 + model.index1(-2.6);
  1461 + test_value(model.index1(), -2.6);
  1462 +
  1463 + // Test breakenergy method
  1464 + model.breakenergy(GEnergy(0.5, "TeV"));
  1465 + test_value(model.breakenergy().TeV(), 0.5);
  1466 +
  1467 + // Test index2 method
  1468 + model.index2(-3.6);
  1469 + test_value(model.index2(), -3.6);
  1470 +
  1471 + // Test operator access
  1472 + const char* strarray[] = {"Prefactor", "Index1", "BreakValue", "Index2"};
  1473 + for (int i = 0; i < 4; ++i) {
  1474 + std::string keyname(strarray[i]);
  1475 + model[keyname].remove_range(); // To allow setting of any value
  1476 + model[keyname].value(2.1);
  1477 + model[keyname].error(1.9);
  1478 + model[keyname].gradient(0.8);
  1479 + test_value(model[keyname].value(), 2.1);
  1480 + test_value(model[keyname].error(), 1.9);
  1481 + test_value(model[keyname].gradient(), 0.8);
  1482 + }
  1483 +
  1484 + // Success if we reached this point
  1485 + test_try_success();
  1486 + }
  1487 + catch (std::exception &e) {
  1488 + test_try_failure(e);
  1489 + }
  1490 +
  1491 + // Exit test
  1492 + return;
  1493 +}
  1494 +
  1495 +
  1496 +/***********************************************************************//**
1409 1497 * @brief Test GModelSpectralLogParabola class
1410 1498 ***************************************************************************/
1411 1499 void TestGModel::test_logparabola(void)
... ... @@ -1812,6 +1900,7 @@ void TestGModel::test_spectral_model(void)
1812 1900 test_xml_model("GModelSpectralPlaw", m_xml_model_point_plaw);
1813 1901 test_xml_model("GModelSpectralPlaw2", m_xml_model_point_plaw2);
1814 1902 test_xml_model("GModelSpectralExpPaw", m_xml_model_point_eplaw);
  1903 + test_xml_model("GModelSpectralBrokenPlaw", m_xml_model_point_bplaw);
1815 1904 test_xml_model("GModelSpectralLogParabola", m_xml_model_point_logparabola);
1816 1905 test_xml_model("GModelSpectralNodes", m_xml_model_point_nodes);
1817 1906 test_xml_model("GModelSpectralFunc", m_xml_model_point_filefct);
... ...
test/test_GModel.hpp
  View file @4960658
... ... @@ -61,6 +61,7 @@ public:
61 61 void test_plaw(void);
62 62 void test_plaw2(void);
63 63 void test_eplaw(void);
  64 + void test_bplaw(void);
64 65 void test_logparabola(void);
65 66 void test_nodes(void);
66 67 void test_filefct(void);
... ... @@ -81,6 +82,7 @@ private:
81 82 std::string m_xml_model_point_plaw;
82 83 std::string m_xml_model_point_plaw2;
83 84 std::string m_xml_model_point_eplaw;
  85 + std::string m_xml_model_point_bplaw;
84 86 std::string m_xml_model_point_logparabola;
85 87 std::string m_xml_model_point_nodes;
86 88 std::string m_xml_model_point_filefct;
... ...