Michael Mayer / gammalib

Commit ea45cec799ff2491bcfff8eed6320f04255ef815

Authored by Alexander Ziegler
1 parent dc225283

Code clean up after review. See change list below. 

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