GammaLib  1.7.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-2018 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())
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  * @exception GException::erange_invalid
573  * Energy range is invalid (emin < emax required).
574  *
575  * Returns Monte Carlo energy by randomly drawing from a smoothly broken
576  * power law.
577  ***************************************************************************/
579  const GEnergy& emax,
580  const GTime& time,
581  GRan& ran) const
582 {
583  // Throw exception if energy range is not valid
584  if (emin >= emax) {
585  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
586  "Minimum energy < maximum energy required.");
587  }
588 
589  // Allocate energy
590  GEnergy energy;
591 
592  // Update Monte Carlo cache
593  update_mc_cache();
594 
595  // Initialse acceptance fraction
596  double acceptance_fraction(0.0);
597 
598  // Use rejection method to draw a random energy. We first draw
599  // analytically from a broken power law, and then compare the power law
600  // at the drawn energy to the curved function. This gives an acceptance
601  // fraction, and we accept the energy only if a uniform random number
602  // is <= the acceptance fraction.
603  do {
604  // Generate an energy from the broken power law distribution
605  energy = m_mc_brokenplaw.mc(emin, emax, time, ran);
606 
607  // Compute acceptance fraction
608  acceptance_fraction = eval(energy) / m_mc_brokenplaw.eval(energy);
609 
610  } while (ran.uniform() > acceptance_fraction);
611 
612  // Return energy
613  return energy;
614 }
615 
616 
617 /***********************************************************************//**
618  * @brief Read model from XML element
619  *
620  * @param[in] xml XML element.
621  *
622  * Reads the spectral information from an XML element.
623  ***************************************************************************/
625 {
626  // Get remaining XML parameters
633 
634  // Read parameters
635  m_norm.read(*prefactor);
636  m_index1.read(*index1);
637  m_index2.read(*index2);
638  m_pivot.read(*pivot);
639  m_breakenergy.read(*breakenergy);
640  m_beta.read(*beta);
641 
642  // Return
643  return;
644 }
645 
646 
647 /***********************************************************************//**
648  * @brief Write model into XML element
649  *
650  * @param[in] xml XML element.
651  *
652  * @exception GException::model_invalid_spectral
653  * Existing XML element is not of the expected type.
654  *
655  * Writes the spectral information into an XML element.
656  ***************************************************************************/
658 {
659  // Set model type
660  if (xml.attribute("type") == "") {
661  xml.attribute("type", type());
662  }
663 
664  // Verify model type
665  if (xml.attribute("type") != type()) {
667  "Spectral model is not of type \""+type()+"\".");
668  }
669 
670  // Get XML parameters
677 
678  // Write parameters
679  m_norm.write(*norm);
680  m_index1.write(*index1);
681  m_pivot.write(*pivot);
682  m_index2.write(*index2);
683  m_breakenergy.write(*breakenergy);
684  m_beta.write(*beta);
685 
686  // Return
687  return;
688 }
689 
690 
691 /***********************************************************************//**
692  * @brief Print smooth broken powerlaw information
693  *
694  * @param[in] chatter Chattiness.
695  * @return String containing model information.
696  ***************************************************************************/
697 std::string GModelSpectralSmoothBrokenPlaw::print(const GChatter& chatter) const
698 {
699  // Initialise result string
700  std::string result;
701 
702  // Continue only if chatter is not silent
703  if (chatter != SILENT) {
704 
705  // Append header
706  result.append("=== GModelSpectralSmoothBrokenPlaw ===");
707 
708  // Append information
709  result.append("\n"+gammalib::parformat("Number of parameters"));
710  result.append(gammalib::str(size()));
711  for (int i = 0; i < size(); ++i) {
712  result.append("\n"+m_pars[i]->print(chatter));
713  }
714 
715  } // endif: chatter was not silent
716 
717  // Return result
718  return result;
719 }
720 
721 
722 /*==========================================================================
723  = =
724  = Private methods =
725  = =
726  ==========================================================================*/
727 
728 /***********************************************************************//**
729  * @brief Initialise class members
730  ***************************************************************************/
732 {
733  // Initialise model type
734  m_type = "SmoothBrokenPowerLaw";
735 
736  // Initialise pre factor
737  m_norm.clear();
738  m_norm.name("Prefactor");
739  m_norm.unit("ph/cm2/s/MeV");
740  m_norm.scale(1.0);
741  m_norm.value(1.0); // default: 1.0
742  m_norm.min(0.0); // min: 0.0
743  m_norm.free();
744  m_norm.gradient(0.0);
745  m_norm.has_grad(true);
746 
747  // Initialise index1
748  m_index1.clear();
749  m_index1.name("Index1");
750  m_index1.scale(1.0);
751  m_index1.value(-2.0); // default: -2.0
752  m_index1.range(-10.0,+10.0); // range: [-10,+10]
753  m_index1.free();
754  m_index1.gradient(0.0);
755  m_index1.has_grad(true);
756 
757  // Initialise index2
758  m_index2.clear();
759  m_index2.name("Index2");
760  m_index2.scale(1.0);
761  m_index2.value(-3.0); // default: -2.0
762  m_index2.range(-10.0,+10.0); // range: [-10,+10]
763  m_index2.free();
764  m_index2.gradient(0.0);
765  m_index2.has_grad(true);
766 
767  // Initialise pivot energy
768  m_pivot.clear();
769  m_pivot.name("PivotEnergy");
770  m_pivot.unit("MeV");
771  m_pivot.scale(1.0e5);
772  m_pivot.value(1.0); // default: 100 GeV
773  m_pivot.fix();
774  m_pivot.gradient(0.0);
775  m_pivot.has_grad(true);
776 
777  // Initialise break energy
779  m_breakenergy.name("BreakEnergy");
780  m_breakenergy.unit("MeV");
781  m_breakenergy.scale(1.0e5);
782  m_breakenergy.value(1.0); // default: 100 GeV
784  m_breakenergy.gradient(0.0);
785  m_breakenergy.has_grad(true);
786 
787  // Initialize beta (break smoothness parameter)
788  m_beta.clear();
789  m_beta.name("BreakSmoothness");
790  m_beta.scale(1.0);
791  m_beta.value(1.0);
792  m_beta.range(0.01,10.0);
793  m_beta.free();
794  m_beta.gradient(0.0);
795  m_beta.has_grad(true);
796 
797  // Set parameter pointer(s)
798  m_pars.clear();
799  m_pars.push_back(&m_norm);
800  m_pars.push_back(&m_index1);
801  m_pars.push_back(&m_pivot);
802  m_pars.push_back(&m_index2);
803  m_pars.push_back(&m_breakenergy);
804  m_pars.push_back(&m_beta);
805 
806  // Initialise eval cache
808  m_last_index1 = 1.0e30;
809  m_last_index2 = 1.0e30;
810  m_last_pivot = 1.0e30;
811  m_last_breakenergy = 1.0e30;
812  m_last_beta = 1.0e30;
813  m_last_epivot_norm = 1.0e30;
814  m_last_ebreak_norm = 1.0e30;
815  m_last_log_epivot_norm = 1.0e30;
816  m_last_log_ebreak_norm = 1.0e30;
817  m_last_epivot_pow = 1.0e30;
818  m_last_ebreak_pow = 1.0e30;
819 
820  // Initialise MC cache
821  m_mc_prefactor = 1.0e30;
822  m_mc_index1 = 1.0e30;
823  m_mc_index2 = 1.0e30;
824  m_mc_pivot = 1.0e30;
825  m_mc_breakenergy = 1.0e30;
827 
828  // Return
829  return;
830 }
831 
832 
833 /***********************************************************************//**
834  * @brief Copy class members
835  *
836  * @param[in] model Smooth broken power law.
837  ***************************************************************************/
839 {
840  // Copy members
841  m_type = model.m_type;
842  m_norm = model.m_norm;
843  m_index1 = model.m_index1;
844  m_index2 = model.m_index2;
845  m_pivot = model.m_pivot;
847  m_beta = model.m_beta;
848 
849  // Set parameter pointer(s)
850  m_pars.clear();
851  m_pars.push_back(&m_norm);
852  m_pars.push_back(&m_index1);
853  m_pars.push_back(&m_pivot);
854  m_pars.push_back(&m_index2);
855  m_pars.push_back(&m_breakenergy);
856  m_pars.push_back(&m_beta);
857 
858  // Copy eval cache
862  m_last_pivot = model.m_last_pivot;
864  m_last_beta = model.m_last_beta;
871 
872  // Copy MC cache
874  m_mc_index1 = model.m_mc_index1;
875  m_mc_index2 = model.m_mc_index2;
876  m_mc_pivot = model.m_mc_pivot;
879 
880  // Return
881  return;
882 }
883 
884 
885 /***********************************************************************//**
886  * @brief Delete class members
887  ***************************************************************************/
889 {
890  // Return
891  return;
892 }
893 
894 
895 /***********************************************************************//**
896  * @brief Update eval precomputation cache
897  *
898  * @param[in] energy Energy.
899  *
900  * Updates the precomputation cache for eval() the method.
901  ***************************************************************************/
903 {
904  // Get parameter values
905  double index1 = m_index1.value();
906  double index2 = m_index2.value();
907  double pivot_eng = pivot().MeV();
908  double break_eng = breakenergy().MeV();
909  double beta = m_beta.value();
910 
911  // If the energy or one of the parameters index1, index2, breakenergy
912  // energy, or beta has changed then recompute the cache
913  if ((m_last_energy != energy) ||
914  (m_last_index1 != index1) ||
915  (m_last_index2 != index2) ||
916  (m_last_pivot != pivot_eng) ||
917  (m_last_breakenergy != break_eng) ||
918  (m_last_beta != beta)) {
919 
920  // Store actual energy and parameter values
921  m_last_energy = energy;
924  m_last_pivot = pivot_eng;
925  m_last_breakenergy = break_eng;
926  m_last_beta = beta;
927 
928  // Compute and store value
929  double eng = energy.MeV();
934 
938 
939  } // endif: recomputation was required
940 
941  // Return
942  return;
943 }
944 
945 
946 /***********************************************************************//**
947  * @brief Update Monte Carlo pre computation cache
948  *
949  * Updates the precomputation cache for Monte Carlo simulations.
950  ***************************************************************************/
952 {
953  // Check if we need to update the cache
954  if (prefactor() != m_mc_prefactor ||
955  index1() != m_mc_index1 ||
956  index2() != m_mc_index2 ||
957  pivot().MeV() != m_mc_pivot ||
958  breakenergy().MeV() != m_mc_breakenergy) {
959 
960  // Set parameters for which Monte Carlo cache was computed
962  m_mc_index1 = index1();
963  m_mc_index2 = index2();
964  m_mc_pivot = pivot().MeV();
966 
967  // Compute prefactor at pivot energy
968  double pre = prefactor() *
969  std::pow(breakenergy().MeV()/pivot().MeV(), index1());
970 
971  // Find out which index is harder. This is important since the smoothly
972  // broken power law follows the hard index below the break energy and
973  // the softer index above the break energy.
976 
977  // Set broken power law for Monte Carlo simulations
979  index1,
980  breakenergy(),
981  index2);
982 
983  } // endif: Update was required
984 
985  // Return
986  return;
987 }
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 gradient factor.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
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:353
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:47
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:54
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
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:184
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
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.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
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:1513
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:1184
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
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.
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:206
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:1332
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:1022
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:1562
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 value factor.
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:413