GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralExpPlaw.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralExpPlaw.cpp - Exponential cut off power law model *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-2016 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GModelSpectralExpPlaw.cpp
23  * @brief Exponential cut off power law spectral class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GException.hpp"
33 #include "GTools.hpp"
34 #include "GRan.hpp"
35 #include "GIntegral.hpp"
38 
39 /* __ Constants __________________________________________________________ */
40 
41 /* __ Globals ____________________________________________________________ */
42 const GModelSpectralExpPlaw g_spectral_eplaw_seed1("ExponentialCutoffPowerLaw",
43  "Prefactor",
44  "Index",
45  "PivotEnergy",
46  "CutoffEnergy");
47 const GModelSpectralRegistry g_spectral_eplaw_registry1(&g_spectral_eplaw_seed1);
48 #if defined(G_LEGACY_XML_FORMAT)
49 const GModelSpectralExpPlaw g_spectral_eplaw_seed2("ExpCutoff",
50  "Prefactor",
51  "Index",
52  "Scale",
53  "Cutoff");
54 const GModelSpectralRegistry g_spectral_eplaw_registry2(&g_spectral_eplaw_seed2);
55 #endif
56 
57 /* __ Method name definitions ____________________________________________ */
58 #define G_MC "GModelSpectralExpPlaw::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
59 #define G_READ "GModelSpectralExpPlaw::read(GXmlElement&)"
60 #define G_WRITE "GModelSpectralExpPlaw::write(GXmlElement&)"
61 
62 /* __ Macros _____________________________________________________________ */
63 
64 /* __ Coding definitions _________________________________________________ */
65 
66 /* __ Debug definitions __________________________________________________ */
67 
68 
69 /*==========================================================================
70  = =
71  = Constructors/destructors =
72  = =
73  ==========================================================================*/
74 
75 /***********************************************************************//**
76  * @brief Void constructor
77  ***************************************************************************/
79 {
80  // Initialise members
81  init_members();
82 
83  // Return
84  return;
85 }
86 
87 
88 /***********************************************************************//**
89  * @brief Model type and parameter name constructor
90  *
91  * @param[in] type Model type.
92  * @param[in] prefactor Name of prefactor parameter.
93  * @param[in] index Name of index parameter.
94  * @param[in] pivot Name of pivot parameter.
95  * @param[in] cutoff Name of cutoff parameter.
96  ***************************************************************************/
98  const std::string& prefactor,
99  const std::string& index,
100  const std::string& pivot,
101  const std::string& cutoff) :
103 {
104  // Initialise members
105  init_members();
106 
107  // Set model type
108  m_type = type;
109 
110  // Set parameter names
111  m_norm.name(prefactor);
112  m_index.name(index);
113  m_pivot.name(pivot);
114  m_ecut.name(cutoff);
115 
116  // Return
117  return;
118 }
119 
120 
121 /***********************************************************************//**
122  * @brief Constructor
123  *
124  * @param[in] prefactor Pre factor normalization (ph/cm2/s/MeV).
125  * @param[in] index Power law index.
126  * @param[in] pivot Pivot energy.
127  * @param[in] cutoff Cut off energy.
128  *
129  * Construct an exponentially cut off power law from
130  * - a prefactor value (in units of ph/cm2/s/MeV)
131  * - a spectral index,
132  * - a pivot energy, and
133  * - a cut off energy.
134  ***************************************************************************/
136  const double& index,
137  const GEnergy& pivot,
138  const GEnergy& cutoff) :
140 {
141  // Initialise members
142  init_members();
143 
144  // Set parameters
145  m_norm.value(prefactor);
146  m_index.value(index);
147  m_pivot.value(pivot.MeV()); // Internally stored in MeV
148  m_ecut.value(cutoff.MeV()); // Internally stored in MeV
149 
150  // Autoscale parameters
151  autoscale();
152 
153  // Return
154  return;
155 }
156 
157 
158 /***********************************************************************//**
159  * @brief XML constructor
160  *
161  * @param[in] xml XML element.
162  *
163  * Constructs an exponentially cut off power law spectral model by extracting
164  * information from an XML element. See the read() method for more
165  * information about the expected structure of the XML element.
166  ***************************************************************************/
169 {
170  // Initialise members
171  init_members();
172 
173  // Read information from XML element
174  read(xml);
175 
176  // Return
177  return;
178 }
179 
180 
181 /***********************************************************************//**
182  * @brief Copy constructor
183  *
184  * @param[in] model Exponentially cut off power law model.
185  ***************************************************************************/
187  GModelSpectral(model)
188 {
189  // Initialise members
190  init_members();
191 
192  // Copy members
193  copy_members(model);
194 
195  // Return
196  return;
197 }
198 
199 
200 /***********************************************************************//**
201  * @brief Destructor
202  ***************************************************************************/
204 {
205  // Free members
206  free_members();
207 
208  // Return
209  return;
210 }
211 
212 
213 /*==========================================================================
214  = =
215  = Operators =
216  = =
217  ==========================================================================*/
218 
219 /***********************************************************************//**
220  * @brief Assignment operator
221  *
222  * @param[in] model Exponentially cut off power law model.
223  * @return Exponentially cut off power law model.
224  ***************************************************************************/
226 {
227  // Execute only if object is not identical
228  if (this != &model) {
229 
230  // Copy base class members
231  this->GModelSpectral::operator=(model);
232 
233  // Free members
234  free_members();
235 
236  // Initialise members
237  init_members();
238 
239  // Copy members
240  copy_members(model);
241 
242  } // endif: object was not identical
243 
244  // Return
245  return *this;
246 }
247 
248 
249 /*==========================================================================
250  = =
251  = Public methods =
252  = =
253  ==========================================================================*/
254 
255 /***********************************************************************//**
256  * @brief Clear exponentially cut off power law model
257  ***************************************************************************/
259 {
260  // Free class members (base and derived classes, derived class first)
261  free_members();
263 
264  // Initialise members
266  init_members();
267 
268  // Return
269  return;
270 }
271 
272 
273 /***********************************************************************//**
274  * @brief Clone exponentially cut off power law model
275  *
276  * @return Pointer to deep copy of exponentially cut off power law model.
277  ***************************************************************************/
279 {
280  // Clone exponentially cut off power law model
281  return new GModelSpectralExpPlaw(*this);
282 }
283 
284 
285 /***********************************************************************//**
286  * @brief Evaluate function
287  *
288  * @param[in] srcEng True photon energy.
289  * @param[in] srcTime True photon arrival time.
290  * @param[in] gradients Compute gradients?
291  * @return Model value (ph/cm2/s/MeV).
292  *
293  * Evaluates
294  *
295  * \f[
296  * S_{\rm E}(E | t) = {\tt m\_norm}
297  * \left( \frac{E}{\tt m\_pivot} \right)^{\tt m\_index}
298  * \exp \left( \frac{-E}{\tt m\_ecut} \right)
299  * \f]
300  *
301  * where
302  * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
303  * - \f${\tt m\_pivot}\f$ is the pivot energy,
304  * - \f${\tt m\_index}\f$ is the spectral index, and
305  * - \f${\tt m\_ecut}\f$ is the cut off energy.
306  *
307  * If the @p gradients flag is true the method will also compute the
308  * partial derivatives of the model with respect to the parameters using
309  *
310  * \f[
311  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} =
312  * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}}
313  * \f]
314  *
315  * \f[
316  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} =
317  * S_{\rm E}(E | t) \, \ln(E/{\tt m_pivot})
318  * \f]
319  *
320  * \f[
321  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_ecut}} =
322  * S_{\rm E}(E | t) \, \left( \frac{E}{{\tt m\_ecut}^2} \right)
323  * \f]
324  *
325  * \f[
326  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_pivot}} =
327  * -S_{\rm E}(E | t) \,
328  * \left( \frac{{\tt m\_index}}{{\tt m\_pivot}} \right)
329  * \f]
330  *
331  * @todo The method expects that pivot!=0 and ecut!=0. Otherwise Inf or NaN
332  * may result. We should add a test that prevents using invalid
333  * values.
334  ***************************************************************************/
336  const GTime& srcTime,
337  const bool& gradients) const
338 {
339  // Update the evaluation cache
340  update_eval_cache(srcEng);
341 
342  // Compute function value
343  double value = m_norm.value() * m_last_power;
344 
345  // Optionally compute gradients
346  if (gradients) {
347 
348  // Compute partial derivatives with respect to the parameter factor
349  // values. The partial derivatives with respect to the parameter
350  // values are obtained by division by the scale factor.
351  double g_norm = (m_norm.is_free())
353  : 0.0;
354  double g_index = (m_index.is_free())
355  ? value * m_index.scale() * std::log(m_last_e_norm)
356  : 0.0;
357  double g_ecut = (m_ecut.is_free())
358  ? value * m_last_e_cut / m_ecut.factor_value()
359  : 0.0;
360  double g_pivot = (m_pivot.is_free())
361  ? -value * m_last_index / m_pivot.factor_value()
362  : 0.0;
363 
364  // Set gradients
365  m_norm.factor_gradient(g_norm);
366  m_index.factor_gradient(g_index);
367  m_ecut.factor_gradient(g_ecut);
368  m_pivot.factor_gradient(g_pivot);
369 
370  } // endif: gradient computation was requested
371 
372  // Compile option: Check for NaN/Inf
373  #if defined(G_NAN_CHECK)
374  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
375  std::cout << "*** ERROR: GModelSpectralExpPlaw::eval";
376  std::cout << "(srcEng=" << srcEng;
377  std::cout << ", srcTime=" << srcTime << "):";
378  std::cout << " NaN/Inf encountered";
379  std::cout << " (value=" << value;
380  std::cout << ", power=" << m_last_power;
381  std::cout << ")" << std::endl;
382  }
383  #endif
384 
385  // Return
386  return value;
387 }
388 
389 
390 /***********************************************************************//**
391  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
392  *
393  * @param[in] emin Minimum photon energy.
394  * @param[in] emax Maximum photon energy.
395  * @return Photon flux (ph/cm2/s).
396  *
397  * Computes
398  *
399  * \f[
400  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
401  * \f]
402  *
403  * where
404  * - [@p emin, @p emax] is an energy interval, and
405  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
406  * The integration is done numerically.
407  ***************************************************************************/
409  const GEnergy& emax) const
410 {
411  // Initialise flux
412  double flux = 0.0;
413 
414  // Compute only if integration range is valid
415  if (emin < emax) {
416 
417  // Setup integration kernel
418  flux_kernel integrand(m_norm.value(), m_index.value(),
419  m_pivot.value(), m_ecut.value());
420  GIntegral integral(&integrand);
421 
422  // Get integration boundaries in MeV
423  double e_min = emin.MeV();
424  double e_max = emax.MeV();
425 
426  // Perform integration
427  flux = integral.romberg(e_min, e_max);
428 
429  } // endif: integration range was valid
430 
431  // Return
432  return flux;
433 }
434 
435 
436 /***********************************************************************//**
437  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
438  *
439  * @param[in] emin Minimum photon energy.
440  * @param[in] emax Maximum photon energy.
441  * @return Energy flux (erg/cm2/s).
442  *
443  * Computes
444  *
445  * \f[
446  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
447  * \f]
448  *
449  * where
450  * - [@p emin, @p emax] is an energy interval, and
451  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
452  * The integration is done numerically.
453  ***************************************************************************/
455  const GEnergy& emax) const
456 {
457  // Initialise flux
458  double eflux = 0.0;
459 
460  // Compute only if integration range is valid
461  if (emin < emax) {
462 
463  // Setup integration kernel
464  eflux_kernel integrand(m_norm.value(), m_index.value(),
465  m_pivot.value(), m_ecut.value());
466  GIntegral integral(&integrand);
467 
468  // Get integration boundaries in MeV
469  double e_min = emin.MeV();
470  double e_max = emax.MeV();
471 
472  // Perform integration
473  eflux = integral.romberg(e_min, e_max);
474 
475  // Convert from MeV/cm2/s to erg/cm2/s
476  eflux *= gammalib::MeV2erg;
477 
478  } // endif: integration range was valid
479 
480  // Return
481  return eflux;
482 }
483 
484 
485 /***********************************************************************//**
486  * @brief Returns MC energy between [emin, emax]
487  *
488  * @param[in] emin Minimum photon energy.
489  * @param[in] emax Maximum photon energy.
490  * @param[in] time True photon arrival time.
491  * @param[in,out] ran Random number generator.
492  * @return Energy.
493  *
494  * @exception GException::erange_invalid
495  * Energy range is invalid (emin < emax required).
496  *
497  * Simulates a random energy in the interval [emin, emax] for an
498  * exponentially cut off power law. The simulation is done using a rejection
499  * method. First, a random energy within [emin, emax] is drawn from an
500  * power law distribution. Then the energy is accepted or rejected based
501  * on an acceptance fraction that is computed from the exponential cut off.
502  ***************************************************************************/
504  const GEnergy& emax,
505  const GTime& time,
506  GRan& ran) const
507 {
508  // Throw an exception if energy range is invalid
509  if (emin >= emax) {
510  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
511  "Minimum energy < maximum energy required.");
512  }
513 
514  // Allocate energy
515  GEnergy energy;
516 
517  // Update cache
518  update_mc_cache(emin, emax);
519 
520  // Initialise energy
521  double eng;
522 
523  // Initialse acceptance fraction
524  double acceptance_fraction;
525  double inv_ecut = 1.0 / m_ecut.value();
526 
527  // Use rejection method to draw a random energy. We first draw
528  // analytically from a power law, and then compare the power law
529  // at the drawn energy to the exponentially cut off function. This
530  // gives an acceptance fraction, and we accept the energy only if
531  // a uniform random number is <= the acceptance fraction.
532  do {
533 
534  // Get uniform random number
535  double u = ran.uniform();
536 
537  // Case A: Index is not -1
538  if (index() != -1.0) {
539  if (u > 0.0) {
541  m_mc_exponent);
542  }
543  else {
544  eng = 0.0;
545  }
546  }
547 
548  // Case B: Index is -1
549  else {
551  }
552 
553  // Compute acceptance fraction
554  acceptance_fraction = std::exp(-eng * inv_ecut);
555 
556  } while (ran.uniform() > acceptance_fraction);
557 
558  // Set energy
559  energy.MeV(eng);
560 
561  // Return energy
562  return energy;
563 }
564 
565 
566 /***********************************************************************//**
567  * @brief Read model from XML element
568  *
569  * @param[in] xml XML element.
570  *
571  * Reads the spectral information from an XML element.
572  ***************************************************************************/
574 {
575  // Get parameter pointers
578  const GXmlElement* ecut = gammalib::xml_get_par(G_READ, xml, m_ecut.name());
580 
581  // Read parameters
582  m_norm.read(*norm);
583  m_index.read(*index);
584  m_ecut.read(*ecut);
585  m_pivot.read(*pivot);
586 
587  // Return
588  return;
589 }
590 
591 
592 /***********************************************************************//**
593  * @brief Write model into XML element
594  *
595  * @param[in] xml XML element.
596  *
597  * @exception GException::model_invalid_spectral
598  * Existing XML element is not of the expected type.
599  *
600  * Writes the spectral information into an XML element.
601  ***************************************************************************/
603 {
604  // Set model type
605  if (xml.attribute("type") == "") {
606  xml.attribute("type", type());
607  }
608 
609  // Verify model type
610  if (xml.attribute("type") != type()) {
612  "Spectral model is not of type \""+type()+"\".");
613  }
614 
615  // Get XML parameters
620 
621  // Write parameters
622  m_norm.write(*norm);
623  m_index.write(*index);
624  m_ecut.write(*ecut);
625  m_pivot.write(*pivot);
626 
627  // Return
628  return;
629 }
630 
631 
632 /***********************************************************************//**
633  * @brief Print model information
634  *
635  * @param[in] chatter Chattiness (defaults to NORMAL).
636  * @return String containing model information.
637  ***************************************************************************/
638 std::string GModelSpectralExpPlaw::print(const GChatter& chatter) const
639 {
640  // Initialise result string
641  std::string result;
642 
643  // Continue only if chatter is not silent
644  if (chatter != SILENT) {
645 
646  // Append header
647  result.append("=== GModelSpectralExpPlaw ===");
648 
649  // Append information
650  result.append("\n"+gammalib::parformat("Number of parameters"));
651  result.append(gammalib::str(size()));
652  for (int i = 0; i < size(); ++i) {
653  result.append("\n"+m_pars[i]->print(chatter));
654  }
655 
656  } // endif: chatter was not silent
657 
658  // Return result
659  return result;
660 }
661 
662 
663 /*==========================================================================
664  = =
665  = Private methods =
666  = =
667  ==========================================================================*/
668 
669 /***********************************************************************//**
670  * @brief Initialise class members
671  ***************************************************************************/
673 {
674  // Initialise model type
675  m_type = "ExponentialCutoffPowerLaw";
676 
677  // Initialise powerlaw normalisation
678  m_norm.clear();
679  m_norm.name("Prefactor");
680  m_norm.unit("ph/cm2/s/MeV");
681  m_norm.scale(1.0);
682  m_norm.value(1.0); // default: 1.0
683  m_norm.min(0.0); // min: 0.0
684  m_norm.free();
685  m_norm.gradient(0.0);
686  m_norm.has_grad(true);
687 
688  // Initialise powerlaw index
689  m_index.clear();
690  m_index.name("Index");
691  m_index.scale(1.0);
692  m_index.value(-2.0); // default: -2.0
693  m_index.range(-10.0,+10.0); // range: [-10,+10]
694  m_index.free();
695  m_index.gradient(0.0);
696  m_index.has_grad(true);
697 
698  // Initialise cut off energy
699  m_ecut.clear();
700  m_ecut.name("CutoffEnergy");
701  m_ecut.unit("MeV");
702  m_ecut.scale(1.0);
703  m_ecut.value(1000.0); // default: 1000.0
704  m_ecut.min(0.1); // min: 0.1
705  m_ecut.free();
706  m_ecut.gradient(0.0);
707  m_ecut.has_grad(true);
708 
709  // Initialise pivot energy
710  m_pivot.clear();
711  m_pivot.name("PivotEnergy");
712  m_pivot.unit("MeV");
713  m_pivot.scale(1.0);
714  m_pivot.value(100.0);
715  m_pivot.fix();
716  m_pivot.gradient(0.0);
717  m_pivot.has_grad(true);
718 
719  // Set parameter pointer(s)
720  m_pars.clear();
721  m_pars.push_back(&m_norm);
722  m_pars.push_back(&m_index);
723  m_pars.push_back(&m_ecut);
724  m_pars.push_back(&m_pivot);
725 
726  // Initialise eval cache
728  m_last_index = 1.0e30;
729  m_last_ecut = 1.0e30;
730  m_last_pivot = 1.0e30;
731  m_last_e_norm = 0.0;
732  m_last_e_cut = 0.0;
733  m_last_power = 0.0;
734 
735  // Initialise MC cache
736  m_mc_emin = 0.0;
737  m_mc_emax = 0.0;
738  m_mc_exponent = 0.0;
739  m_mc_pow_emin = 0.0;
740  m_mc_pow_ewidth = 0.0;
741 
742  // Return
743  return;
744 }
745 
746 
747 /***********************************************************************//**
748  * @brief Copy class members
749  *
750  * @param[in] model Exponential cut off power law model.
751  ***************************************************************************/
753 {
754  // Copy members
755  m_type = model.m_type;
756  m_norm = model.m_norm;
757  m_index = model.m_index;
758  m_ecut = model.m_ecut;
759  m_pivot = model.m_pivot;
760 
761  // Set parameter pointer(s)
762  m_pars.clear();
763  m_pars.push_back(&m_norm);
764  m_pars.push_back(&m_index);
765  m_pars.push_back(&m_ecut);
766  m_pars.push_back(&m_pivot);
767 
768  // Copy eval cache
770  m_last_index = model.m_last_index;
771  m_last_ecut = model.m_last_ecut;
772  m_last_pivot = model.m_last_pivot;
774  m_last_e_cut = model.m_last_e_cut;
775  m_last_power = model.m_last_power;
776 
777  // Copy MC cache
778  m_mc_emin = model.m_mc_emin;
779  m_mc_emax = model.m_mc_emax;
783 
784  // Return
785  return;
786 }
787 
788 
789 /***********************************************************************//**
790  * @brief Delete class members
791  ***************************************************************************/
793 {
794  // Return
795  return;
796 }
797 
798 
799 /***********************************************************************//**
800  * @brief Update eval precomputation cache
801  *
802  * @param[in] energy Energy.
803  *
804  * Updates the precomputation cache for eval() method.
805  ***************************************************************************/
807 {
808  // Get parameter values (takes 3 multiplications which are difficult
809  // to avoid)
810  double index = m_index.value();
811  double ecut = m_ecut.value();
812  double pivot = m_pivot.value();
813 
814  // If the energy or one of the parameters index, cut-off or pivot
815  // energy has changed then recompute the cache
816  if ((m_last_energy != energy) ||
817  (m_last_index != index) ||
818  (m_last_ecut != ecut) ||
819  (m_last_pivot != pivot)) {
820 
821  // Store actual energy and parameter values
822  m_last_energy = energy;
824  m_last_ecut = ecut;
826 
827  // Compute and store value
828  double eng = energy.MeV();
829  m_last_e_norm = eng / m_last_pivot;
830  m_last_e_cut = eng / m_last_ecut;
833  } // endif: recomputation was required
834 
835  // Return
836  return;
837 }
838 
839 
840 /***********************************************************************//**
841  * @brief Update Monte Carlo pre computation cache
842  *
843  * @param[in] emin Minimum photon energy.
844  * @param[in] emax Maximum photon energy.
845  *
846  * Updates the precomputation cache for Monte Carlo simulations.
847  ***************************************************************************/
849  const GEnergy& emax) const
850 
851 {
852  // Case A: Index is not -1
853  if (index() != -1.0) {
854 
855  // Change in energy boundaries?
856  if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) {
857  m_mc_emin = emin.MeV();
858  m_mc_emax = emax.MeV();
859  m_mc_exponent = index() + 1.0;
862  }
863 
864  }
865 
866  // Case B: Index is -1
867  else {
868 
869  // Change in energy boundaries?
870  if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) {
871  m_mc_emin = emin.MeV();
872  m_mc_emax = emax.MeV();
873  m_mc_exponent = 0.0;
876  }
877 
878  }
879 
880  // Return
881  return;
882 }
883 
884 
885 /***********************************************************************//**
886  * @brief Kernel for photon flux integration
887  *
888  * @param[in] energy Energy (MeV).
889  ***************************************************************************/
890 double GModelSpectralExpPlaw::flux_kernel::eval(const double& energy)
891 {
892  // Evaluate function value
893  double e_norm = energy * m_inv_pivot;
894  double e_cut = energy * m_inv_ecut;
895  double power = std::pow(e_norm, m_index) * std::exp(-e_cut);
896  double value = m_norm * power;
897 
898  // Return value
899  return value;
900 }
901 
902 
903 /***********************************************************************//**
904  * @brief Kernel for energy flux integration
905  *
906  * @param[in] energy Energy (MeV).
907  ***************************************************************************/
908 double GModelSpectralExpPlaw::eflux_kernel::eval(const double& energy)
909 {
910  // Evaluate function value
911  double e_norm = energy * m_inv_pivot;
912  double e_cut = energy * m_inv_ecut;
913  double power = std::pow(e_norm, m_index) * std::exp(-e_cut);
914  double value = m_norm * power * energy;
915 
916  // Return value
917  return value;
918 }
const double & factor_gradient(void) const
Return parameter gradient factor.
double m_last_e_norm
Last E/Epivot value.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
#define G_MC
double m_mc_emax
Maximum energy.
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
Abstract spectral model base class.
double m_last_index
Last index parameter.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
virtual std::string type(void) const
Return model type.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
virtual void write(GXmlElement &xml) const
Write model into XML element.
void copy_members(const GModelSpectralExpPlaw &model)
Copy class members.
std::vector< GModelPar * > m_pars
Parameter pointers.
double m_mc_exponent
Exponent (index+1)
GEnergy pivot(void) const
Return pivot energy.
XML element node class.
Definition: GXmlElement.hpp:47
Spectral model registry class definition.
double eval(const double &eng)
Kernel for photon flux integration.
Random number generator class.
Definition: GRan.hpp:44
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
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.
GModelPar m_ecut
Exponential cut off energy.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
const GModelSpectralExpPlaw g_spectral_eplaw_seed1("ExponentialCutoffPowerLaw","Prefactor","Index","PivotEnergy","CutoffEnergy")
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
virtual ~GModelSpectralExpPlaw(void)
Destructor.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
GModelPar m_norm
Normalization factor.
virtual void read(const GXmlElement &xml)
Read model from XML element.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns MC energy between [emin, emax].
const GXmlAttribute * attribute(const int &index) const
Return attribute.
virtual void clear(void)
Clear exponentially cut off power law model.
void free(void)
Free a parameter.
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
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1184
double m_last_ecut
Last energy cut-off parameter.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
double m_mc_emin
Minimum energy.
void clear(void)
Clear parameter.
GModelSpectralExpPlaw(void)
Void constructor.
double index(void) const
Return power law index.
GChatter
Definition: GTypemaps.hpp:33
GEnergy m_last_energy
Last energy value.
std::string m_type
Model type.
double m_last_pivot
Last pivot parameter.
Interface definition for the spectral model registry class.
void autoscale(void)
Autoscale parameters.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print model information.
double m_last_power
Last power value.
void init_members(void)
Initialise class members.
double m_mc_pow_emin
Power of minimum energy.
const double MeV2erg
Definition: GTools.hpp:45
void update_mc_cache(const GEnergy &emin, const GEnergy &emax) const
Update Monte Carlo pre computation cache.
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.
void update_eval_cache(const GEnergy &energy) const
Update eval precomputation cache.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
const std::string & unit(void) const
Return parameter unit.
Exception handler interface definition.
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
Exponential cut off power law spectral class.
#define G_READ
double eval(const double &eng)
Kernel for energy flux integration.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
GModelPar m_pivot
Pivot energy.
virtual GModelSpectralExpPlaw & operator=(const GModelSpectralExpPlaw &model)
Assignment operator.
void init_members(void)
Initialise class members.
Exponential cut off power law spectral class interface definition.
virtual GModelSpectralExpPlaw * clone(void) const
Clone exponentially cut off power law model.
#define G_WRITE
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
GModelPar m_index
Spectral index.
Integration class interface definition.
double m_last_e_cut
Last E/Ecut value.
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 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
double m_mc_pow_ewidth
Power of energy width.