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