GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralSmoothBrokenPlaw.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralSmoothBrokenPlaw.hpp *
3  * Smoothly broken power law spectrum class *
4  * ----------------------------------------------------------------------- *
5  * copyright (C) 2017-2018 by Josh Cardenzana *
6  * ----------------------------------------------------------------------- *
7  * *
8  * This program is free software: you can redistribute it and/or modify *
9  * it under the terms of the GNU General Public License as published by *
10  * the Free Software Foundation, either version 3 of the License, or *
11  * (at your option) any later version. *
12  * *
13  * This program is distributed in the hope that it will be useful, *
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16  * GNU General Public License for more details. *
17  * *
18  * You should have received a copy of the GNU General Public License *
19  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
20  * *
21  ***************************************************************************/
22 /**
23  * @file GModelSpectralSmoothBrokenPlaw.hpp
24  * @brief Smoothly broken power law spectrum class definition
25  * @author Josh Cardenzana
26  */
27 
28 #ifndef GMODELSPECTRALSMOOTHBROKENPLAW_HPP
29 #define GMODELSPECTRALSMOOTHBROKENPLAW_HPP
30 
31 /* __ Includes ___________________________________________________________ */
32 #include <string>
33 #include "GModelSpectral.hpp"
34 #include "GFunction.hpp"
35 #include "GModelPar.hpp"
36 #include "GEnergy.hpp"
38 
39 /* __ Forward declarations _______________________________________________ */
40 class GRan;
41 class GTime;
42 class GXmlElement;
43 
44 
45 /***********************************************************************//**
46  * @class GModelSpectralSmoothBrokenPlaw
47  *
48  * @brief Smoothly broken power law spectral model class
49  *
50  * This class implements a smoothly broken power law spectrum. The model is
51  * defined by
52  *
53  * \f[
54  * S_{\rm E}(E | t) = k_0 \left( \frac{E}{E_0} \right)^{\gamma_1}
55  * \left[ 1 + \left( \frac{E}{E_b} \right)^{\frac{\gamma_1 - \gamma_2}{\beta}}
56  * \right]^{-\beta}
57  * \f]
58  *
59  * where
60  * - \f$k_0\f$ is the normalization or prefactor,
61  * - \f$\gamma_1\f$ is the spectral index before the break,
62  * - \f$\gamma_2\f$ is the spectral index after the break,
63  * - \f$E_0\f$ is the pivot energy,
64  * - \f$E_b\f$ is the break energy,
65  * - \f$\beta\f$ defines the break smoothness (larger = smoother transition)
66  *
67  * Notes:
68  * - The pivot energy should be set far away from the expected break energy
69  * value.
70  * - When the two indices are close together, the beta parameter becomes
71  * poorly constrained. Since the beta parameter also scales the indices, this
72  * can cause very large errors in the estimates of the various spectral
73  * parameters. In this case, consider fixing beta.
74  ***************************************************************************/
76 
77 public:
78  // Constructors and destructors
80  GModelSpectralSmoothBrokenPlaw(const std::string& type,
81  const std::string& prefactor,
82  const std::string& index1,
83  const std::string& pivot,
84  const std::string& index2,
85  const std::string& breakenergy,
86  const std::string& beta);
87  GModelSpectralSmoothBrokenPlaw(const double& prefactor,
88  const double& index1,
89  const GEnergy& pivot,
90  const double& index2,
91  const GEnergy& breakenergy,
92  const double& beta);
93  explicit GModelSpectralSmoothBrokenPlaw(const GXmlElement& xml);
95  virtual ~GModelSpectralSmoothBrokenPlaw(void);
96 
97  // Operators
99 
100  // Implemented pure virtual methods
101  virtual void clear(void);
102  virtual GModelSpectralSmoothBrokenPlaw* clone(void) const;
103  virtual std::string classname(void) const;
104  virtual std::string type(void) const;
105  virtual double eval(const GEnergy& srcEng,
106  const GTime& srcTime = GTime(),
107  const bool& gradients = false) const;
108  virtual double flux(const GEnergy& emin,
109  const GEnergy& emax) const;
110  virtual double eflux(const GEnergy& emin,
111  const GEnergy& emax) const;
112  virtual GEnergy mc(const GEnergy& emin,
113  const GEnergy& emax,
114  const GTime& time,
115  GRan& ran) const;
116  virtual void read(const GXmlElement& xml);
117  virtual void write(GXmlElement& xml) const;
118  virtual std::string print(const GChatter& chatter = NORMAL) const;
119 
120  // Other methods
121  void type(const std::string& type);
122  double prefactor(void) const;
123  double index1(void) const;
124  double index2(void) const;
125  GEnergy pivot(void) const;
126  GEnergy breakenergy(void) const;
127  double beta(void) const;
128  void prefactor(const double& prefactor);
129  void index1(const double& index1);
130  void index2(const double& index2);
131  void pivot(const GEnergy& pivot);
132  void breakenergy(const GEnergy& breakenergy);
133  void beta(const double& beta);
134 
135 protected:
136  // Protected methods
137  void init_members(void);
138  void copy_members(const GModelSpectralSmoothBrokenPlaw& model);
139  void free_members(void);
140  void update_eval_cache(const GEnergy& energy) const;
141  void update_mc_cache(void) const;
142 
143  // Class to determine the integral photon flux
144  class flux_kern : public GFunction {
145  public:
146  // Constructor
147  flux_kern(const double& prefactor,
148  const double& index1,
149  const GEnergy& pivot,
150  const double& index2,
151  const GEnergy& breakenergy,
152  const double& beta) :
153  m_prefactor(prefactor),
154  m_index1(index1),
155  m_index2(index2),
156  m_pivot(pivot),
157  m_breakenergy(breakenergy),
158  m_beta(beta) {};
159 
160  // Evaluation method
161  double eval(const double& energy) {
162  double epivot = energy / m_pivot.MeV();
163  double ebreak = energy / m_breakenergy.MeV();
164  return m_prefactor * std::pow(epivot, m_index1) *
165  std::pow(1.0 + std::pow(ebreak,(m_index1-m_index2)/m_beta),-m_beta);
166  }
167 
168  // Members
169  protected:
170  double m_prefactor; //!< Normalization
171  double m_index1; //!< Spectral index1
172  double m_index2; //!< Spectral index2
173  GEnergy m_pivot; //!< Pivot energy
174  GEnergy m_breakenergy; //!< Break energy
175  double m_beta; //!< Break smoothness parameter
176  };
177 
178  // Class to determine the integral energy flux, derivation of flux_kern
179  class eflux_kern : public flux_kern {
180  public:
181  // Constructor
182  eflux_kern(const double& prefactor,
183  const double& index1,
184  const GEnergy& pivot,
185  const double& index2,
186  const GEnergy& breakenergy,
187  const double& beta):
188  flux_kern(prefactor, index1, pivot, index2, breakenergy, beta)
189  {};
190 
191  // Evaluation method
192  double eval(const double& energy) {
193  return energy * flux_kern::eval(energy);
194  }
195  };
196 
197  // Protected members
198  std::string m_type; //!< Model type
199  GModelPar m_norm; //!< Normalization factor
200  GModelPar m_index1; //!< Spectral index1
201  GModelPar m_index2; //!< Spectral index2
202  GModelPar m_pivot; //!< Pivot energy
203  GModelPar m_breakenergy; //!< Energy of spectral break
204  GModelPar m_beta; //!< Break smoothness
205 
206  // Cached members used for pre-computations
207  mutable GEnergy m_last_energy; //!< Last energy value
208  mutable double m_last_index1; //!< Last index1 parameter
209  mutable double m_last_index2; //!< Last index2 parameter
210  mutable double m_last_pivot; //!< Last pivot parameter
211  mutable double m_last_breakenergy; //!< Last breakenergy parameter
212  mutable double m_last_beta; //!< Last beta parameter
213  mutable double m_last_epivot_norm; //!< Last E/Epivot value
214  mutable double m_last_ebreak_norm; //!< Last E/Ebreakenergy value
215  mutable double m_last_log_epivot_norm; //!< Last ln(E/Epivot) value
216  mutable double m_last_log_ebreak_norm; //!< Last ln(E/Ebreakenergy) value
217  mutable double m_last_epivot_pow; //!< Last pow(E/Epivot,index1) value
218  mutable double m_last_ebreak_pow; //!< Last pow(E/Ebreakenergy,(index1-index2)/beta)
219 
220  // Cached members for Monte-Carlo simulations
221  mutable double m_mc_prefactor; //!< Last pre factor
222  mutable double m_mc_index1; //!< Last first index
223  mutable double m_mc_index2; //!< Last second index
224  mutable double m_mc_pivot; //!< Last pivot energy
225  mutable double m_mc_breakenergy; //!< Last break energy
226  mutable GModelSpectralBrokenPlaw m_mc_brokenplaw; //!< Broken power plaw
227 };
228 
229 
230 /***********************************************************************//**
231  * @brief Return class name
232  *
233  * @return String containing the class name ("GModelSpectralSmoothBrokenPlaw").
234  ***************************************************************************/
235 inline
237 {
238  return ("GModelSpectralSmoothBrokenPlaw");
239 }
240 
241 
242 /***********************************************************************//**
243  * @brief Return model type
244  *
245  * @return "SmoothBrokenPowerLaw".
246  *
247  * Returns the type of the spectral smoothly broken power law model.
248  ***************************************************************************/
249 inline
251 {
252  return (m_type);
253 }
254 
255 
256 /***********************************************************************//**
257  * @brief Set model type
258  *
259  * @param[in] type Model type.
260  *
261  * Set the type of the spectral smoothly broken power law model.
262  ***************************************************************************/
263 inline
264 void GModelSpectralSmoothBrokenPlaw::type(const std::string& type)
265 {
266  m_type = type;
267  return;
268 }
269 
270 
271 /***********************************************************************//**
272  * @brief Return pre factor
273  *
274  * @return Pre factor (ph/cm2/s/MeV).
275  *
276  * Returns the pre factor.
277  ***************************************************************************/
278 inline
280 {
281  return (m_norm.value());
282 }
283 
284 
285 /***********************************************************************//**
286  * @brief Set pre factor
287  *
288  * @param[in] prefactor Pre factor (ph/cm2/s/MeV).
289  *
290  * Sets the pre factor.
291  ***************************************************************************/
292 inline
293 void GModelSpectralSmoothBrokenPlaw::prefactor(const double& prefactor)
294 {
295  m_norm.value(prefactor);
296  return;
297 }
298 
299 
300 /***********************************************************************//**
301  * @brief Return smoothly broken power law index1
302  *
303  * @return Power law index1.
304  *
305  * Returns the power law index1.
306  ***************************************************************************/
307 inline
309 {
310  return (m_index1.value());
311 }
312 
313 
314 /***********************************************************************//**
315  * @brief Set smoothly broken power law index1
316  *
317  * @param[in] index1 Power law index1.
318  *
319  * Sets the power law index1.
320  ***************************************************************************/
321 inline
322 void GModelSpectralSmoothBrokenPlaw::index1(const double& index1)
323 {
324  m_index1.value(index1);
325  return;
326 }
327 
328 
329 /***********************************************************************//**
330  * @brief Return smoothly broken power law index2
331  *
332  * @return Power law index2.
333  *
334  * Returns the power law index2.
335  ***************************************************************************/
336 inline
338 {
339  return (m_index2.value());
340 }
341 
342 
343 /***********************************************************************//**
344  * @brief Set smoothly broken power law index2
345  *
346  * @param[in] index2 Power law index2.
347  *
348  * Sets the power law index2.
349  ***************************************************************************/
350 inline
351 void GModelSpectralSmoothBrokenPlaw::index2(const double& index2)
352 {
353  m_index2.value(index2);
354  return;
355 }
356 
357 
358 /***********************************************************************//**
359 * @brief Return pivot energy
360 *
361 * @return Smoothly broken power law pivot energy.
362 *
363 * Returns the smoothly broken power law scale energy.
364 ***************************************************************************/
365 inline
367 {
368  GEnergy energy;
369  energy.MeV(m_pivot.value());
370  return energy;
371 }
372 
373 
374 /***********************************************************************//**
375 * @brief Set pivot energy
376 *
377 * @param[in] pivot Smoothly broken power law pivot energy.
378 *
379 * Sets the smoothly broken power law pivot energy.
380 ***************************************************************************/
381 inline
383 {
384  m_pivot.value(pivot.MeV());
385  return;
386 }
387 
388 
389 /***********************************************************************//**
390  * @brief Return breakenergy energy
391  *
392  * @return breakenergy energy.
393  *
394  * Returns the breakenergy energy.
395  ***************************************************************************/
396 inline
398 {
399  GEnergy energy;
400  energy.MeV(m_breakenergy.value());
401  return energy;
402 }
403 
404 
405 /***********************************************************************//**
406  * @brief Set breakenergy energy
407  *
408  * @param[in] breakenergy break energy.
409  *
410  * Sets the "breakenergy" energy.
411  ***************************************************************************/
412 inline
414 {
415  m_breakenergy.value(breakenergy.MeV());
416  return;
417 }
418 
419 
420 /***********************************************************************//**
421 * @brief Returns break smoothness parameter beta
422 *
423 * @return break smoothness parameter.
424 *
425 * Returns the break smoothness parameter 'beta'.
426 ***************************************************************************/
427 inline
429 {
430  return (m_beta.value());
431 }
432 
433 
434 /***********************************************************************//**
435 * @brief Set break smoothness
436 *
437 * @param[in] beta break smoothness parameter.
438 *
439 * Sets the beta break smoothness parameter.
440 ***************************************************************************/
441 inline
442 void GModelSpectralSmoothBrokenPlaw::beta(const double& beta)
443 {
444  m_beta.value(beta);
445  return;
446 }
447 #endif /* GMODELSPECTRALSMOOTHBROKENPLAW_HPP */
double m_last_ebreak_pow
Last pow(E/Ebreakenergy,(index1-index2)/beta)
void init_members(void)
Initialise class members.
Energy value class definition.
Abstract spectral model base class.
GModelPar m_breakenergy
Energy of spectral break.
virtual ~GModelSpectralSmoothBrokenPlaw(void)
Destructor.
XML element node class.
Definition: GXmlElement.hpp:47
double prefactor(void) const
Return pre factor.
double m_last_ebreak_norm
Last E/Ebreakenergy value.
Random number generator class.
Definition: GRan.hpp:44
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Smoothly broken power law spectral model class.
Time class.
Definition: GTime.hpp:54
void free_members(void)
Delete class members.
void copy_members(const GModelSpectralSmoothBrokenPlaw &model)
Copy class members.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns Monte Carlo energy between [emin, emax].
virtual std::string classname(void) const
Return class name.
Model parameter class interface definition.
Model parameter class.
Definition: GModelPar.hpp:87
double m_last_log_epivot_norm
Last ln(E/Epivot) value.
double index1(void) const
Return smoothly broken power law index1.
Single parameter function abstract base class definition.
double m_last_epivot_pow
Last pow(E/Epivot,index1) value.
GModelSpectralSmoothBrokenPlaw(void)
Void constructor.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print smooth broken powerlaw information.
Broken power law spectrum class definition.
GChatter
Definition: GTypemaps.hpp:33
virtual GModelSpectralSmoothBrokenPlaw * clone(void) const
Clone smoothly broken power law model.
virtual void read(const GXmlElement &xml)
Read model from XML element.
virtual GModelSpectralSmoothBrokenPlaw & operator=(const GModelSpectralSmoothBrokenPlaw &model)
Assignment operator.
Broken power law spectral model class.
GEnergy pivot(void) const
Return pivot energy.
virtual void clear(void)
Clear smoothly broken power law model.
Abstract spectral model base class interface definition.
double index2(void) const
Return smoothly broken power law index2.
GEnergy breakenergy(void) const
Return breakenergy energy.
virtual std::string type(void) const
Return model type.
double value(void) const
Return parameter value.
Single parameter function abstract base class.
Definition: GFunction.hpp:44
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
double m_last_log_ebreak_norm
Last ln(E/Ebreakenergy) value.
double beta(void) const
Returns break smoothness parameter beta.
flux_kern(const double &prefactor, const double &index1, const GEnergy &pivot, const double &index2, const GEnergy &breakenergy, const double &beta)
double m_last_breakenergy
Last breakenergy parameter.
GModelSpectralBrokenPlaw m_mc_brokenplaw
Broken power plaw.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
virtual void write(GXmlElement &xml) const
Write model into XML element.
eflux_kern(const double &prefactor, const double &index1, const GEnergy &pivot, const double &index2, const GEnergy &breakenergy, const double &beta)
void update_eval_cache(const GEnergy &energy) const
Update eval precomputation cache.
void update_mc_cache(void) const
Update Monte Carlo pre computation cache.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48