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