GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralSmoothBrokenPlaw.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralSmoothBrokenPlaw.cpp *
3  * Smoothly broken power law spectrum class *
4  * ----------------------------------------------------------------------- *
5  * copyright (C) 2017-2021 by Joshua 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.cpp
24  * @brief Smoothly broken power law spectrum class implementation
25  * @author Joshua Cardenzana
26  */
27 
28 /* __ Includes ___________________________________________________________ */
29 #ifdef HAVE_CONFIG_H
30 #include <config.h>
31 #endif
32 #include <cmath>
33 #include "GException.hpp"
34 #include "GIntegral.hpp"
35 #include "GTools.hpp"
36 #include "GRan.hpp"
39 
40 /* __ Constants __________________________________________________________ */
41 
42 /* __ Globals ____________________________________________________________ */
43 const GModelSpectralSmoothBrokenPlaw g_spectral_sblaw_seed1("SmoothBrokenPowerLaw",
44  "Prefactor",
45  "Index1",
46  "PivotEnergy",
47  "Index2",
48  "BreakEnergy",
49  "BreakSmoothness");
50 const GModelSpectralRegistry g_spectral_sblaw_registry1(&g_spectral_sblaw_seed1);
51 const GModelSpectralSmoothBrokenPlaw g_spectral_sblaw_seed2("SmoothBrokenPowerLaw",
52  "Prefactor",
53  "Index1",
54  "Scale",
55  "Index2",
56  "BreakValue",
57  "Beta");
58 const GModelSpectralRegistry g_spectral_sblaw_registry2(&g_spectral_sblaw_seed2);
59 
60 /* __ Method name definitions ____________________________________________ */
61 #define G_MC "GModelSpectralSmoothBrokenPlaw::mc(GEnergy&, GEnergy&, GTime&,"\
62  " GRan&)"
63 #define G_READ "GModelSpectralSmoothBrokenPlaw::read(GXmlElement&)"
64 #define G_WRITE "GModelSpectralSmoothBrokenPlaw::write(GXmlElement&)"
65 
66 /* __ Macros _____________________________________________________________ */
67 
68 /* __ Coding definitions _________________________________________________ */
69 
70 /* __ Debug definitions __________________________________________________ */
71 
72 
73 /*==========================================================================
74  = =
75  = Constructors/destructors =
76  = =
77  ==========================================================================*/
78 
79 /***********************************************************************//**
80  * @brief Void constructor
81  ***************************************************************************/
84 {
85  // Initialise private members for clean destruction
86  init_members();
87 
88  // Return
89  return;
90 }
91 
92 
93 /***********************************************************************//**
94  * @brief Model type and parameter name constructor
95  *
96  * @param[in] type Model type.
97  * @param[in] prefactor Name of prefactor parameter.
98  * @param[in] index1 Name of index1 parameter.
99  * @param[in] pivot Name of pivot parameter.
100  * @param[in] index2 Name of index2 parameter.
101  * @param[in] breakenergy Name of breakenergy parameter.
102  * @param[in] beta Name of beta parameter.
103  ***************************************************************************/
105  const std::string& type,
106  const std::string& prefactor,
107  const std::string& index1,
108  const std::string& pivot,
109  const std::string& index2,
110  const std::string& breakenergy,
111  const std::string& beta) :
113 {
114  // Initialise members
115  init_members();
116 
117  // Set model type
118  m_type = type;
119 
120  // Set parameter names
121  m_norm.name(prefactor);
122  m_index1.name(index1);
123  m_pivot.name(pivot);
124  m_index2.name(index2);
125  m_breakenergy.name(breakenergy);
126  m_beta.name(beta);
127 
128  // Return
129  return;
130 }
131 
132 
133 /***********************************************************************//**
134  * @brief Parameter constructor
135  *
136  * @param[in] prefactor Smoothly broken power law pre factor (ph/cm2/s/MeV).
137  * @param[in] index1 Smoothly broken power law index1.
138  * @param[in] pivot Smoothly broken power law pivot energy
139  * @param[in] index2 Smoothly broken power law index1.
140  * @param[in] breakenergy Break energy.
141  * @param[in] beta Break smoothness parameter
142  *
143  * Constructs a smoothly broken power law using the model parameters
144  * power law @p prefactor (ph/cm2/s/MeV),
145  * spectral @p index1,
146  * @p pivot energy,
147  * spectral @p index2,
148  * @p breakenergy of spectral break, and
149  * smoothness parameter @p beta.
150  ***************************************************************************/
152  const double& prefactor,
153  const double& index1,
154  const GEnergy& pivot,
155  const double& index2,
156  const GEnergy& breakenergy,
157  const double& beta) :
159 {
160  // Initialise members
161  init_members();
162 
163  // Set parameters
164  m_norm.value(prefactor);
165  m_index1.value(index1);
166  m_pivot.value(pivot.MeV()); // Internally stored in MeV
167  m_index2.value(index2);
168  m_breakenergy.value(breakenergy.MeV()); // Internally stored in MeV
169  m_beta.value(beta);
170 
171  // Perform autoscaling of parameter
172  autoscale();
173 
174  // Return
175  return;
176 }
177 
178 
179 /***********************************************************************//**
180  * @brief XML constructor
181  *
182  * @param[in] xml XML element containing position information.
183  *
184  * Constructs a smoothly broken power law spectral model by extracting
185  * information from an XML element. See the read() method for more
186  * information about the expected structure of the XML element.
187  ***************************************************************************/
190 {
191  // Initialise members
192  init_members();
193 
194  // Read information from XML element
195  read(xml);
196 
197  // Return
198  return;
199 }
200 
201 
202 /***********************************************************************//**
203  * @brief Copy constructor
204  *
205  * @param[in] model Smoothly broken power law model.
206  ***************************************************************************/
208  const GModelSpectralSmoothBrokenPlaw& model) :
209  GModelSpectral(model)
210 {
211  // Initialise members
212  init_members();
213 
214  // Copy members
215  copy_members(model);
216 
217  // Return
218  return;
219 }
220 
221 
222 /***********************************************************************//**
223  * @brief Destructor
224  ***************************************************************************/
226 {
227  // Free members
228  free_members();
229 
230  // Return
231  return;
232 }
233 
234 
235 /*==========================================================================
236  = =
237  = Operators =
238  = =
239  ==========================================================================*/
240 
241 /***********************************************************************//**
242  * @brief Assignment operator
243  *
244  * @param[in] model Smoothly broken power law model.
245  * @return Smoothly broken power law model.
246  ***************************************************************************/
248  const GModelSpectralSmoothBrokenPlaw& model)
249 {
250  // Execute only if object is not identical
251  if (this != &model) {
252 
253  // Copy base class members
254  this->GModelSpectral::operator=(model);
255 
256  // Free members
257  free_members();
258 
259  // Initialise members
260  init_members();
261 
262  // Copy members
263  copy_members(model);
264 
265  } // endif: object was not identical
266 
267  // Return
268  return *this;
269 }
270 
271 
272 /*==========================================================================
273  = =
274  = Public methods =
275  = =
276  ==========================================================================*/
277 
278 /***********************************************************************//**
279  * @brief Clear smoothly broken power law model
280  ***************************************************************************/
282 {
283  // Free class members (base and derived classes, derived class first)
284  free_members();
286 
287  // Initialise members
289  init_members();
290 
291  // Return
292  return;
293 }
294 
295 
296 /***********************************************************************//**
297  * @brief Clone smoothly broken power law model
298  *
299  * @return Pointer to deep copy of smoothly broken power law model.
300  ***************************************************************************/
302 {
303  // Clone spectral broken power law model
304  return new GModelSpectralSmoothBrokenPlaw(*this);
305 }
306 
307 
308 /***********************************************************************//**
309  * @brief Evaluate function
310  *
311  * @param[in] srcEng True photon energy.
312  * @param[in] srcTime True photon arrival time.
313  * @param[in] gradients Compute gradients?
314  * @return Model value (ph/cm2/s/MeV).
315  *
316  * Evaluates
317  *
318  * \f[
319  * S_{\rm E}(E | t) = k_0 \left( \frac{E}{E_0} \right)^{\gamma_1}
320  * \left[ 1 + \left( \frac{E}{E_b} \right)^{\frac{\gamma_1 - \gamma_2}{\beta}}
321  * \right]^{-\beta}
322  * \f]
323  *
324  * where:
325  * - \f$k_0\f$ is the normalization or prefactor,
326  * - \f$\gamma_1\f$ is the spectral index before the break,
327  * - \f$\gamma_2\f$ is the spectral index after the break,
328  * - \f$E_0\f$ is the pivot energy,
329  * - \f$E_b\f$ is the break energy,
330  * - \f$\beta\f$ is the break smoothness.
331  *
332  * If the @p gradients flag is true the method will also compute the
333  * partial derivatives of the model with respect to the parameters using
334  *
335  * \f[
336  * \frac{\delta S_{\rm E}(E | t)}{\delta k_0} =
337  * \frac{S_{\rm E}(E | t)}{k_0}
338  * \f]
339  *
340  * \f[
341  * \frac{\delta S_{\rm E}(E | t)}{\delta \gamma_1} =
342  * S_{\rm E}(E | t) \left[ \ln\left( \frac{E}{E_0} \right) -
343  * \frac{\left(\frac{E}{E_b}\right)^{\frac{\gamma_1 - \gamma_2}{\beta}}
344  * \ln\left( \frac{E}{E_b} \right)}
345  * {\left(\frac{E}{E_b}\right)^{\frac{\gamma_1 - \gamma_2}{\beta}} + 1}\right]
346  * \f]
347  *
348  * \f[
349  * \frac{\delta S_{\rm E}(E | t)}{\delta \gamma_2} =
350  * S_{\rm E}(E | t)
351  * \frac{\left(\frac{E}{E_b}\right)^{\frac{\gamma_1 - \gamma_2}{\beta}}
352  * \ln\left( \frac{E}{E_b} \right)}
353  * {\left(\frac{E}{E_b}\right)^{\frac{\gamma_1 - \gamma_2}{\beta}} + 1}
354  * \f]
355  *
356  * \f[
357  * \frac{\delta S_{\rm E}(E | t)}{\delta E_0} =
358  * S_{\rm E}(E | t) \frac{\gamma_1}{E_0}
359  * \f]
360  *
361  * \f[
362  * \frac{\delta S_{\rm E}(E | t)}{\delta E_b} =
363  * S_{\rm E}(E | t) \frac{\left(\gamma_1 - \gamma_2 \right)
364  * \left( \frac{E}{E_b} \right)^{\frac{\gamma_1 - \gamma_2}{\beta}}}
365  {E_b \left(1 + \left(\frac{E}{E_b}\right)^{\frac{\gamma_1 - \gamma_2}{\beta}} \right)}
366  * \f]
367  *
368  * \f[
369  * \frac{\delta S_{\rm E}(E | t)}{\delta \beta} = S_{\rm E}(E | t)
370  * \left[
371  * \frac{ \left( \frac{E}{E_b} \right)^{\frac{\gamma_1 - \gamma_2}{\beta}}
372  * \ln \left( \left( \frac{E}{E_b}\right)^{ \frac{\gamma_1 - \gamma_2}{\beta}} \right)}
373  * {1 + \left(\frac{E}{E_b}\right)^{\frac{\gamma_1 - \gamma_2}{\beta}}}
374  * - \ln \left( 1 + \left(\frac{E}{E_b}\right)^{\frac{\gamma_1 - \gamma_2}{\beta}} \right)
375  * \right]
376  * \f]
377  *
378  * @todo The method expects that energy!=0. Otherwise Inf or NaN may result.
379  ***************************************************************************/
381  const GTime& srcTime,
382  const bool& gradients) const
383 {
384  // Update the evaluation cache
385  update_eval_cache(srcEng);
386 
387  // Compute function value
388  double value = m_norm.value() * m_last_epivot_pow *
390 
391  // Optionally compute gradients
392  if (gradients) {
393 
394  // Compute normalisation gradient
395  double g_norm = (m_norm.is_free())
398  : 0.0;
399 
400  // Compute index1 and index2 value gradients
401  double g_index1 = (m_index1.is_free())
402  ? value * m_index1.scale() * (m_last_log_epivot_norm -
404  (m_last_ebreak_pow + 1.0)))
405  : 0.0;
406  double g_index2 = (m_index2.is_free())
407  ? value * m_index2.scale() * m_last_log_ebreak_norm *
409  : 0.0;
410 
411  // Compute pivot and break energy value gradients
412  double g_pivot = (m_pivot.is_free() && m_pivot.factor_value() != 0.0)
413  ? -value * m_last_index1 / m_pivot.factor_value()
414  : 0.0;
415  double g_break = (m_breakenergy.is_free())
418  : 0.0;
419 
420  // Compute beta gradient
421  double g_beta = (m_beta.is_free())
422  ? value * m_beta.scale() *
425  : 0.0;
426 
427  // Store the gradient values
428  m_norm.factor_gradient(g_norm);
429  m_index1.factor_gradient(g_index1);
430  m_index2.factor_gradient(g_index2);
431  m_pivot.factor_gradient(g_pivot);
433  m_beta.factor_gradient(g_beta);
434 
435  } // endif: gradient computation was requested
436 
437  // Compile option: Check for NaN/Inf
438  #if defined(G_NAN_CHECK)
439  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
440  std::cout << "*** ERROR: GModelSpectralSmoothBrokenPlaw::eval";
441  std::cout << "(srcEng=" << srcEng;
442  std::cout << ", srcTime=" << srcTime << "):";
443  std::cout << " NaN/Inf encountered";
444  std::cout << " (value=" << value;
445  std::cout << ", m_norm=" << m_norm.value();
446  std::cout << ", m_index1=" << m_index1.value();
447  std::cout << ", m_index2=" << m_index2.value();
448  std::cout << ", m_pivot=" << m_pivot.value();
449  std::cout << ", m_breakenergy=" << m_breakenergy.value();
450  std::cout << ", m_beta=" << m_beta.value();
451  std::cout << ", m_last_epivot_norm=" << m_last_epivot_norm;
452  std::cout << ", m_last_ebreak_norm=" << m_last_ebreak_norm;
453  std::cout << ", m_last_log_epivot_norm=" << m_last_log_epivot_norm;
454  std::cout << ", m_last_log_ebreak_norm=" << m_last_log_ebreak_norm;
455  std::cout << ", m_last_epivot_pow=" << m_last_epivot_pow;
456  std::cout << ", m_last_ebreak_pow=" << m_last_ebreak_pow;
457  std::cout << ")" << std::endl;
458  }
459  #endif
460 
461  // Return
462  return value;
463 }
464 
465 
466 /***********************************************************************//**
467  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
468  *
469  * @param[in] emin Minimum photon energy.
470  * @param[in] emax Maximum photon energy.
471  * @return Photon flux (ph/cm2/s).
472  *
473  * Computes
474  *
475  * \f[
476  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
477  * \f]
478  *
479  * where
480  * - [@p emin, @p emax] is an energy interval, and
481  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
482  * The integration is done numerically.
483  ***************************************************************************/
485  const GEnergy& emax) const
486 {
487  // Initialise flux
488  double flux = 0.0;
489 
490  // Compute only if integration range is valid
491  if (emin < emax) {
492 
493  // Initialise function to integrate
494  flux_kern kernel(prefactor(), index1(), pivot(),
495  index2(), breakenergy(), beta());
496 
497  // Initialise integral class with function
498  GIntegral integral(&kernel);
499 
500  // Set integration precision
501  integral.eps(1.0e-8);
502 
503  // Calculate integral between emin and emax
504  flux = integral.romberg(emin.MeV(), emax.MeV());
505 
506  } // endif: integration range was valid
507 
508  // Return flux
509  return flux;
510 }
511 
512 
513 /***********************************************************************//**
514  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
515  *
516  * @param[in] emin Minimum photon energy.
517  * @param[in] emax Maximum photon energy.
518  * @return Energy flux (erg/cm2/s).
519  *
520  * Computes
521  *
522  * \f[
523  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
524  * \f]
525  *
526  * where
527  * - [@p emin, @p emax] is an energy interval, and
528  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
529  * The integration is done numerically.
530  ***************************************************************************/
532  const GEnergy& emax) const
533 {
534  // Initialise flux
535  double eflux = 0.0;
536 
537  // Compute only if integration range is valid
538  if (emin < emax) {
539 
540  // Initialise function to integrate
541  eflux_kern kernel(prefactor(), index1(), pivot(),
542  index2(), breakenergy(), beta());
543 
544  // Initialise integral class with function
545  GIntegral integral(&kernel);
546 
547  // Set integration precision
548  integral.eps(1.0e-8);
549 
550  // Calculate integral between emin and emax
551  eflux = integral.romberg(emin.MeV(), emax.MeV());
552 
553  // Convert from MeV/cm2/s to erg/cm2/s
554  eflux *= gammalib::MeV2erg;
555 
556  } // endif: integration range was valid
557 
558  // Return flux
559  return eflux;
560 }
561 
562 
563 /***********************************************************************//**
564  * @brief Returns Monte Carlo energy between [emin, emax]
565  *
566  * @param[in] emin Minimum photon energy.
567  * @param[in] emax Maximum photon energy.
568  * @param[in] time True photon arrival time.
569  * @param[in,out] ran Random number generator.
570  * @return Energy.
571  *
572  * Returns Monte Carlo energy by randomly drawing from a smoothly broken
573  * power law.
574  ***************************************************************************/
576  const GEnergy& emax,
577  const GTime& time,
578  GRan& ran) const
579 {
580  // Check energy interval
582 
583  // Allocate energy
584  GEnergy energy;
585 
586  // Update Monte Carlo cache
587  update_mc_cache();
588 
589  // Initialse acceptance fraction
590  double acceptance_fraction(0.0);
591 
592  // Use rejection method to draw a random energy. We first draw
593  // analytically from a broken power law, and then compare the power law
594  // at the drawn energy to the curved function. This gives an acceptance
595  // fraction, and we accept the energy only if a uniform random number
596  // is <= the acceptance fraction.
597  do {
598  // Generate an energy from the broken power law distribution
599  energy = m_mc_brokenplaw.mc(emin, emax, time, ran);
600 
601  // Compute acceptance fraction
602  acceptance_fraction = eval(energy) / m_mc_brokenplaw.eval(energy);
603 
604  } while (ran.uniform() > acceptance_fraction);
605 
606  // Return energy
607  return energy;
608 }
609 
610 
611 /***********************************************************************//**
612  * @brief Read model from XML element
613  *
614  * @param[in] xml XML element.
615  *
616  * Reads the spectral information from an XML element.
617  ***************************************************************************/
619 {
620  // Verify number of model parameters
622 
623  // Get remaining XML parameters
630 
631  // Read parameters
632  m_norm.read(*prefactor);
633  m_index1.read(*index1);
634  m_index2.read(*index2);
635  m_pivot.read(*pivot);
636  m_breakenergy.read(*breakenergy);
637  m_beta.read(*beta);
638 
639  // Return
640  return;
641 }
642 
643 
644 /***********************************************************************//**
645  * @brief Write model into XML element
646  *
647  * @param[in] xml XML element.
648  *
649  * Writes the spectral information into an XML element.
650  ***************************************************************************/
652 {
653  // Verify model type
655 
656  // Get XML parameters
663 
664  // Write parameters
665  m_norm.write(*norm);
666  m_index1.write(*index1);
667  m_pivot.write(*pivot);
668  m_index2.write(*index2);
669  m_breakenergy.write(*breakenergy);
670  m_beta.write(*beta);
671 
672  // Return
673  return;
674 }
675 
676 
677 /***********************************************************************//**
678  * @brief Print smooth broken powerlaw information
679  *
680  * @param[in] chatter Chattiness.
681  * @return String containing model information.
682  ***************************************************************************/
683 std::string GModelSpectralSmoothBrokenPlaw::print(const GChatter& chatter) const
684 {
685  // Initialise result string
686  std::string result;
687 
688  // Continue only if chatter is not silent
689  if (chatter != SILENT) {
690 
691  // Append header
692  result.append("=== GModelSpectralSmoothBrokenPlaw ===");
693 
694  // Append information
695  result.append("\n"+gammalib::parformat("Number of parameters"));
696  result.append(gammalib::str(size()));
697  for (int i = 0; i < size(); ++i) {
698  result.append("\n"+m_pars[i]->print(chatter));
699  }
700 
701  } // endif: chatter was not silent
702 
703  // Return result
704  return result;
705 }
706 
707 
708 /*==========================================================================
709  = =
710  = Private methods =
711  = =
712  ==========================================================================*/
713 
714 /***********************************************************************//**
715  * @brief Initialise class members
716  ***************************************************************************/
718 {
719  // Initialise model type
720  m_type = "SmoothBrokenPowerLaw";
721 
722  // Initialise pre factor
723  m_norm.clear();
724  m_norm.name("Prefactor");
725  m_norm.unit("ph/cm2/s/MeV");
726  m_norm.scale(1.0);
727  m_norm.value(1.0); // default: 1.0
728  m_norm.min(0.0); // min: 0.0
729  m_norm.free();
730  m_norm.gradient(0.0);
731  m_norm.has_grad(true);
732 
733  // Initialise index1
734  m_index1.clear();
735  m_index1.name("Index1");
736  m_index1.scale(1.0);
737  m_index1.value(-2.0); // default: -2.0
738  m_index1.range(-10.0,+10.0); // range: [-10,+10]
739  m_index1.free();
740  m_index1.gradient(0.0);
741  m_index1.has_grad(true);
742 
743  // Initialise index2
744  m_index2.clear();
745  m_index2.name("Index2");
746  m_index2.scale(1.0);
747  m_index2.value(-3.0); // default: -2.0
748  m_index2.range(-10.0,+10.0); // range: [-10,+10]
749  m_index2.free();
750  m_index2.gradient(0.0);
751  m_index2.has_grad(true);
752 
753  // Initialise pivot energy
754  m_pivot.clear();
755  m_pivot.name("PivotEnergy");
756  m_pivot.unit("MeV");
757  m_pivot.scale(1.0e5);
758  m_pivot.value(1.0); // default: 100 GeV
759  m_pivot.fix();
760  m_pivot.gradient(0.0);
761  m_pivot.has_grad(true);
762 
763  // Initialise break energy
765  m_breakenergy.name("BreakEnergy");
766  m_breakenergy.unit("MeV");
767  m_breakenergy.scale(1.0e5);
768  m_breakenergy.value(1.0); // default: 100 GeV
770  m_breakenergy.gradient(0.0);
771  m_breakenergy.has_grad(true);
772 
773  // Initialize beta (break smoothness parameter)
774  m_beta.clear();
775  m_beta.name("BreakSmoothness");
776  m_beta.scale(1.0);
777  m_beta.value(1.0);
778  m_beta.range(0.01,10.0);
779  m_beta.free();
780  m_beta.gradient(0.0);
781  m_beta.has_grad(true);
782 
783  // Set parameter pointer(s)
784  m_pars.clear();
785  m_pars.push_back(&m_norm);
786  m_pars.push_back(&m_index1);
787  m_pars.push_back(&m_pivot);
788  m_pars.push_back(&m_index2);
789  m_pars.push_back(&m_breakenergy);
790  m_pars.push_back(&m_beta);
791 
792  // Initialise eval cache
794  m_last_index1 = 1.0e30;
795  m_last_index2 = 1.0e30;
796  m_last_pivot = 1.0e30;
797  m_last_breakenergy = 1.0e30;
798  m_last_beta = 1.0e30;
799  m_last_epivot_norm = 1.0e30;
800  m_last_ebreak_norm = 1.0e30;
801  m_last_log_epivot_norm = 1.0e30;
802  m_last_log_ebreak_norm = 1.0e30;
803  m_last_epivot_pow = 1.0e30;
804  m_last_ebreak_pow = 1.0e30;
805 
806  // Initialise MC cache
807  m_mc_prefactor = 1.0e30;
808  m_mc_index1 = 1.0e30;
809  m_mc_index2 = 1.0e30;
810  m_mc_pivot = 1.0e30;
811  m_mc_breakenergy = 1.0e30;
813 
814  // Return
815  return;
816 }
817 
818 
819 /***********************************************************************//**
820  * @brief Copy class members
821  *
822  * @param[in] model Smooth broken power law.
823  ***************************************************************************/
825 {
826  // Copy members
827  m_type = model.m_type;
828  m_norm = model.m_norm;
829  m_index1 = model.m_index1;
830  m_index2 = model.m_index2;
831  m_pivot = model.m_pivot;
833  m_beta = model.m_beta;
834 
835  // Set parameter pointer(s)
836  m_pars.clear();
837  m_pars.push_back(&m_norm);
838  m_pars.push_back(&m_index1);
839  m_pars.push_back(&m_pivot);
840  m_pars.push_back(&m_index2);
841  m_pars.push_back(&m_breakenergy);
842  m_pars.push_back(&m_beta);
843 
844  // Copy eval cache
848  m_last_pivot = model.m_last_pivot;
850  m_last_beta = model.m_last_beta;
857 
858  // Copy MC cache
860  m_mc_index1 = model.m_mc_index1;
861  m_mc_index2 = model.m_mc_index2;
862  m_mc_pivot = model.m_mc_pivot;
865 
866  // Return
867  return;
868 }
869 
870 
871 /***********************************************************************//**
872  * @brief Delete class members
873  ***************************************************************************/
875 {
876  // Return
877  return;
878 }
879 
880 
881 /***********************************************************************//**
882  * @brief Update eval precomputation cache
883  *
884  * @param[in] energy Energy.
885  *
886  * Updates the precomputation cache for eval() the method.
887  ***************************************************************************/
889 {
890  // Get parameter values
891  double index1 = m_index1.value();
892  double index2 = m_index2.value();
893  double pivot_eng = pivot().MeV();
894  double break_eng = breakenergy().MeV();
895  double beta = m_beta.value();
896 
897  // If the energy or one of the parameters index1, index2, breakenergy
898  // energy, or beta has changed then recompute the cache
899  if ((m_last_energy != energy) ||
900  (m_last_index1 != index1) ||
901  (m_last_index2 != index2) ||
902  (m_last_pivot != pivot_eng) ||
903  (m_last_breakenergy != break_eng) ||
904  (m_last_beta != beta)) {
905 
906  // Store actual energy and parameter values
907  m_last_energy = energy;
910  m_last_pivot = pivot_eng;
911  m_last_breakenergy = break_eng;
912  m_last_beta = beta;
913 
914  // Compute and store value
915  double eng = energy.MeV();
920 
924 
925  } // endif: recomputation was required
926 
927  // Return
928  return;
929 }
930 
931 
932 /***********************************************************************//**
933  * @brief Update Monte Carlo pre computation cache
934  *
935  * Updates the precomputation cache for Monte Carlo simulations.
936  ***************************************************************************/
938 {
939  // Check if we need to update the cache
940  if (prefactor() != m_mc_prefactor ||
941  index1() != m_mc_index1 ||
942  index2() != m_mc_index2 ||
943  pivot().MeV() != m_mc_pivot ||
944  breakenergy().MeV() != m_mc_breakenergy) {
945 
946  // Set parameters for which Monte Carlo cache was computed
948  m_mc_index1 = index1();
949  m_mc_index2 = index2();
950  m_mc_pivot = pivot().MeV();
952 
953  // Compute prefactor at pivot energy
954  double pre = prefactor() *
955  std::pow(breakenergy().MeV()/pivot().MeV(), index1());
956 
957  // Find out which index is harder. This is important since the smoothly
958  // broken power law follows the hard index below the break energy and
959  // the softer index above the break energy.
962 
963  // Set broken power law for Monte Carlo simulations
965  index1,
966  breakenergy(),
967  index2);
968 
969  } // endif: Update was required
970 
971  // Return
972  return;
973 }
Smoothly broken power law spectrum class definition.
const GModelSpectralSmoothBrokenPlaw g_spectral_sblaw_seed1("SmoothBrokenPowerLaw","Prefactor","Index1","PivotEnergy","Index2","BreakEnergy","BreakSmoothness")
double m_last_ebreak_pow
Last pow(E/Ebreakenergy,(index1-index2)/beta)
void init_members(void)
Initialise class members.
const GModelSpectralSmoothBrokenPlaw g_spectral_sblaw_seed2("SmoothBrokenPowerLaw","Prefactor","Index1","Scale","Index2","BreakValue","Beta")
const double & factor_gradient(void) const
Return parameter factor gradient.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:382
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
Abstract spectral model base class.
double gradient(void) const
Return parameter gradient.
GModelPar m_breakenergy
Energy of spectral break.
int size(void) const
Return number of parameters.
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
virtual ~GModelSpectralSmoothBrokenPlaw(void)
Destructor.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:48
double prefactor(void) const
Return pre factor.
Spectral model registry class definition.
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:55
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
void free_members(void)
Delete class members.
double min(void) const
Return parameter minimum boundary.
bool is_free(void) const
Signal if parameter is free.
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].
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
double m_last_log_epivot_norm
Last ln(E/Epivot) value.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
double index1(void) const
Return smoothly broken power law index1.
void free(void)
Free a parameter.
double m_last_epivot_pow
Last pow(E/Epivot,index1) value.
void fix(void)
Fix a parameter.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1637
GModelSpectralSmoothBrokenPlaw(void)
Void constructor.
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 print(const GChatter &chatter=NORMAL) const
Print smooth broken powerlaw information.
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void check_energy_interval(const std::string &origin, const GEnergy &emin, const GEnergy &emax)
Checks energy interval.
Definition: GException.cpp:364
void clear(void)
Clear parameter.
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.
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
Interface definition for the spectral model registry class.
virtual void clear(void)
Clear smoothly broken power law model.
void autoscale(void)
Autoscale parameters.
double index2(void) const
Return smoothly broken power law index2.
GEnergy breakenergy(void) const
Return breakenergy energy.
void eps(const double &eps)
Set relative precision.
Definition: GIntegral.hpp:214
void init_members(void)
Initialise class members.
virtual std::string type(void) const
Return model type.
const double MeV2erg
Definition: GTools.hpp:45
void range(const double &min, const double &max)
Set minimum and maximum parameter boundaries.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
double value(void) const
Return parameter value.
const std::string & unit(void) const
Return parameter unit.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
double m_last_log_ebreak_norm
Last ln(E/Ebreakenergy) value.
Exception handler interface definition.
double beta(void) const
Returns break smoothness parameter 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
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
virtual void write(GXmlElement &xml) const
Write model into XML element.
Integration class interface definition.
void update_eval_cache(const GEnergy &energy) const
Update eval precomputation cache.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1689
void free_members(void)
Delete class members.
void update_mc_cache(void) const
Update Monte Carlo pre computation cache.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
const double & factor_value(void) const
Return parameter factor value.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819