GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralLogParabola.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralLogParabola.cpp - Log parabola spectral model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-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 GModelSpectralLogParabola.cpp
23  * @brief Log parabola spectral model class definition
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 ____________________________________________________________ */
43  "Prefactor",
44  "Index",
45  "PivotEnergy",
46  "Curvature");
47 const GModelSpectralRegistry g_spectral_logp_registry1(&g_spectral_logp_seed1);
48 #if defined(G_LEGACY_XML_FORMAT)
49 const GModelSpectralLogParabola g_spectral_logp_seed2("LogParabola",
50  "norm",
51  "alpha",
52  "Eb",
53  "beta");
54 const GModelSpectralLogParabola g_spectral_logp_seed3("LogParabola",
55  "Prefactor",
56  "Index",
57  "Scale",
58  "Curvature");
59 const GModelSpectralRegistry g_spectral_logp_registry2(&g_spectral_logp_seed2);
60 const GModelSpectralRegistry g_spectral_logp_registry3(&g_spectral_logp_seed3);
61 #endif
62 
63 /* __ Method name definitions ____________________________________________ */
64 #define G_FLUX "GModelSpectralLogParabola::flux(GEnergy&, GEnergy&)"
65 #define G_EFLUX "GModelSpectralLogParabola::eflux(GEnergy&, GEnergy&)"
66 #define G_MC "GModelSpectralLogParabola::mc(GEnergy&, GEnergy&, GTime&,"\
67  " GRan&)"
68 #define G_READ "GModelSpectralLogParabola::read(GXmlElement&)"
69 #define G_WRITE "GModelSpectralLogParabola::write(GXmlElement&)"
70 
71 /* __ Macros _____________________________________________________________ */
72 
73 /* __ Coding definitions _________________________________________________ */
74 
75 /* __ Debug definitions __________________________________________________ */
76 
77 
78 /*==========================================================================
79  = =
80  = Constructors/destructors =
81  = =
82  ==========================================================================*/
83 
84 /***********************************************************************//**
85  * @brief Void constructor
86  ***************************************************************************/
88 {
89  // Initialise private members
90  init_members();
91 
92  // Return
93  return;
94 }
95 
96 
97 /***********************************************************************//**
98  * @brief Model type and parameter name constructor
99  *
100  * @param[in] type Model type.
101  * @param[in] prefactor Name of prefactor parameter.
102  * @param[in] index Name of index parameter.
103  * @param[in] pivot Name of pivot parameter.
104  * @param[in] curvature Name of curvature parameter.
105  ***************************************************************************/
107  const std::string& prefactor,
108  const std::string& index,
109  const std::string& pivot,
110  const std::string& curvature) :
112 {
113  // Initialise members
114  init_members();
115 
116  // Set model type
117  m_type = type;
118 
119  // Set parameter names
120  m_norm.name(prefactor);
121  m_index.name(index);
122  m_pivot.name(pivot);
123  m_curvature.name(curvature);
124 
125  // Return
126  return;
127 }
128 
129 
130 /***********************************************************************//**
131  * @brief Constructor
132  *
133  * @param[in] prefactor Power law pre factor.
134  * @param[in] index Power law index.
135  * @param[in] pivot Pivot energy.
136  * @param[in] curvature Curvature.
137  *
138  * Construct a LogParabola model from
139  * - a pre factor (in ph/cm2/s/MeV),
140  * - a spectral index,
141  * - a pivot energy, and
142  * - a curvature.
143  ***************************************************************************/
145  const double& index,
146  const GEnergy& pivot,
147  const double& curvature) :
149 {
150  // Initialise members
151  init_members();
152 
153  // Set parameters
154  m_norm.value(prefactor);
155  m_index.value(index);
156  m_pivot.value(pivot.MeV());
157  m_curvature.value(curvature);
158 
159  // Autoscale parameters
160  autoscale();
161 
162  // Return
163  return;
164 }
165 
166 
167 /***********************************************************************//**
168  * @brief XML constructor
169  *
170  * @param[in] xml XML element.
171  *
172  * Constructs log parabola spectral model by extracting information from an
173  * XML element. See the read() method for more information about the expected
174  * structure of the XML element.
175  ***************************************************************************/
178 {
179  // Initialise members
180  init_members();
181 
182  // Read information from XML element
183  read(xml);
184 
185  // Return
186  return;
187 }
188 
189 
190 /***********************************************************************//**
191  * @brief Copy constructor
192  *
193  * @param[in] model LogParabola model.
194  ***************************************************************************/
196  GModelSpectral(model)
197 {
198  // Initialise members
199  init_members();
200 
201  // Copy members
202  copy_members(model);
203 
204  // Return
205  return;
206 }
207 
208 
209 /***********************************************************************//**
210  * @brief Destructor
211  ***************************************************************************/
213 {
214  // Free members
215  free_members();
216 
217  // Return
218  return;
219 }
220 
221 
222 /*==========================================================================
223  = =
224  = Operators =
225  = =
226  ==========================================================================*/
227 
228 /***********************************************************************//**
229  * @brief Assignment operator
230  *
231  * @param[in] model LogParabola model.
232  * @return LogParabola model.
233  ***************************************************************************/
235 {
236  // Execute only if object is not identical
237  if (this != &model) {
238 
239  // Copy base class members
240  this->GModelSpectral::operator=(model);
241 
242  // Free members
243  free_members();
244 
245  // Initialise members
246  init_members();
247 
248  // Copy members
249  copy_members(model);
250 
251  } // endif: object was not identical
252 
253  // Return
254  return *this;
255 }
256 
257 
258 /*==========================================================================
259  = =
260  = Public methods =
261  = =
262  ==========================================================================*/
263 
264 /***********************************************************************//**
265  * @brief Clear log parabola model
266  ***************************************************************************/
268 {
269  // Free class members (base and derived classes, derived class first)
270  free_members();
272 
273  // Initialise members
275  init_members();
276 
277  // Return
278  return;
279 }
280 
281 
282 /***********************************************************************//**
283  * @brief Clone log parabola model
284  *
285  * @return Pointer to deep copy of log parabola spectrum.
286  ***************************************************************************/
288 {
289  // Clone log parabola model
290  return new GModelSpectralLogParabola(*this);
291 }
292 
293 
294 /***********************************************************************//**
295  * @brief Evaluate model
296  *
297  * @param[in] srcEng True photon energy.
298  * @param[in] srcTime True photon arrival time.
299  * @param[in] gradients Compute gradients?
300  * @return Model value (ph/cm2/s/MeV).
301  *
302  * Computes
303  *
304  * \f[
305  * S_{\rm E}(E | t) = {\tt m\_norm}
306  * \left( \frac{E}{\tt m\_pivot} \right)^{{\tt m\_index} +
307  * {\tt m\_curvature} \, \ln \frac{E}{\tt m\_pivot}}
308  * \f]
309  *
310  * where
311  * - \f${\tt m\_norm}\f$ is the normalization or prefactor,
312  * - \f${\tt m\_index}\f$ is the spectral index,
313  * - \f${\tt m\_curvature}\f$ is the spectral curvature, and
314  * - \f${\tt m\_pivot}\f$ is the pivot energy.
315  *
316  * If the @p gradients flag is true the method will also compute the
317  * partial derivatives of the model with respect to the parameters using
318  *
319  * \f[
320  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_norm}} =
321  * \frac{S_{\rm E}(E | t)}{{\tt m\_norm}}
322  * \f]
323  *
324  * \f[
325  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_index}} =
326  * S_{\rm E}(E | t) \, \ln(E/{\tt m_pivot})
327  * \f]
328  *
329  * \f[
330  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_curvature}} =
331  * S_{\rm E}(E | t) \, (\ln(E/{\tt m_pivot})^2)
332  * \f]
333  *
334  * \f[
335  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_pivot}} =
336  * -S_{\rm E}(E | t) \,
337  * \left( \frac{2 {\tt m\_curvature} \ln(E/{\tt m_pivot}) + {\tt m\_index}}
338  * {{\tt m\_pivot}} \right)
339  * \f]
340  *
341  * @todo The method expects that energy!=0. Otherwise Inf or NaN may result.
342  ***************************************************************************/
344  const GTime& srcTime,
345  const bool& gradients) const
346 {
347  // Update the evaluation cache
348  update_eval_cache(srcEng);
349 
350  // Compute function value
351  double value = m_norm.value() * m_last_power;
352 
353  // Optionally compute gradients
354  if (gradients) {
355 
356  // Compute partial derivatives of the parameter values
357  double g_norm = (m_norm.is_free())
358  ? m_norm.scale() * m_last_power : 0.0;
359  double g_index = (m_index.is_free())
360  ? value * m_index.scale() * m_last_log_e_norm : 0.0;
361  double g_curvature = (m_curvature.is_free())
362  ? value * m_curvature.scale() * m_last_log_e_norm *
363  m_last_log_e_norm : 0.0;
364  double g_pivot = (m_pivot.is_free())
365  ? -value * (m_last_exponent + m_curvature.value() *
367 
368  // Set gradients
369  m_norm.factor_gradient(g_norm);
370  m_index.factor_gradient(g_index);
371  m_curvature.factor_gradient(g_curvature);
372  m_pivot.factor_gradient(g_pivot);
373 
374  } // endif: gradient computation was requested
375 
376  // Compile option: Check for NaN/Inf
377  #if defined(G_NAN_CHECK)
378  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
379  std::cout << "*** ERROR: GModelSpectralLogParabola::eval";
380  std::cout << "(srcEng=" << srcEng;
381  std::cout << ", srcTime=" << srcTime << "):";
382  std::cout << " NaN/Inf encountered";
383  std::cout << " (value=" << value;
384  std::cout << ", power=" << m_last_power;
385  std::cout << ")" << std::endl;
386  }
387  #endif
388 
389  // Return
390  return value;
391 }
392 
393 
394 /***********************************************************************//**
395  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
396  *
397  * @param[in] emin Minimum photon energy.
398  * @param[in] emax Maximum photon energy.
399  * @return Photon flux (ph/cm2/s).
400  *
401  * Computes
402  *
403  * \f[
404  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
405  * \f]
406  *
407  * where
408  * - [@p emin, @p emax] is an energy interval, and
409  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
410  * The integration is done numerically.
411  ***************************************************************************/
413  const GEnergy& emax) const
414 {
415  // Initialise flux
416  double flux = 0.0;
417 
418  // Compute only if integration range is valid
419  if (emin < emax) {
420 
421  // Initialise function to integrate
422  flux_kern kernel(prefactor(), index(), curvature(), pivot());
423 
424  // Initialise integral class with function
425  GIntegral integral(&kernel);
426 
427  // Set integration precision
428  integral.eps(1.0e-8);
429 
430  // Calculate integral between emin and emax
431  flux = integral.romberg(emin.MeV(), emax.MeV());
432 
433  } // endif: integration range was valid
434 
435  // Return flux
436  return flux;
437 }
438 
439 
440 /***********************************************************************//**
441  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
442  *
443  * @param[in] emin Minimum photon energy.
444  * @param[in] emax Maximum photon energy.
445  * @return Energy flux (erg/cm2/s).
446  *
447  * Computes
448  *
449  * \f[
450  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
451  * \f]
452  *
453  * where
454  * - [@p emin, @p emax] is an energy interval, and
455  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
456  * The integration is done numerically.
457  ***************************************************************************/
459  const GEnergy& emax) const
460 {
461  // Initialise flux
462  double eflux = 0.0;
463 
464  // Compute only if integration range is valid
465  if (emin < emax) {
466 
467  // Initialise function to integrate
468  eflux_kern kernel(prefactor(), index(), curvature(), pivot());
469 
470  // Initialise integral class with function
471  GIntegral integral(&kernel);
472 
473  // Set integration precision
474  integral.eps(1.0e-8);
475 
476  // Calculate integral between emin and emax
477  eflux = integral.romberg(emin.MeV(), emax.MeV());
478 
479  // Convert from MeV/cm2/s to erg/cm2/s
480  eflux *= gammalib::MeV2erg;
481 
482  } // endif: integration range was valid
483 
484  // Return flux
485  return eflux;
486 }
487 
488 
489 /***********************************************************************//**
490  * @brief Returns Monte Carlo energy between [emin, emax]
491  *
492  * @param[in] emin Minimum photon energy.
493  * @param[in] emax Maximum photon energy.
494  * @param[in] time True photon arrival time.
495  * @param[in,out] ran Random number generator.
496  * @return Energy.
497  *
498  * @exception GException::erange_invalid
499  * Energy range is invalid (emin < emax required).
500  *
501  * Returns Monte Carlo energy by randomly drawing from the spectral model.
502  ***************************************************************************/
504  const GEnergy& emax,
505  const GTime& time,
506  GRan& ran) const
507 {
508  if (emin >= emax) {
509  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
510  "Minimum energy < maximum energy required.");
511  }
512 
513  // Allocate energy
514  GEnergy energy;
515 
516  // Update cache
517  update_mc_cache(emin, emax, time);
518 
519  // Initialise energy
520  double eng;
521 
522  // Initialse acceptance fraction
523  double acceptance_fraction;
524 
525  // Use rejection method to draw a random energy. We first draw
526  // analytically from a power law, and then compare the power law
527  // at the drawn energy to the curved function. This
528  // gives an acceptance fraction, and we accept the energy only if
529  // a uniform random number is <= the acceptance fraction.
530  do {
531 
532  // Get uniform random number
533  double u = ran.uniform();
534 
535  // Case A: Corresponding mc Plaw-Index is not -1
536  if (m_mc_exponent != 0.0) {
537  if (u > 0.0) {
539  m_mc_exponent);
540  }
541  else {
542  eng = 0.0;
543  }
544  }
545 
546  // Case B: Corresponding mc Plaw-Index is -1
547  else {
549  }
550 
551  // Compute powerlaw at given energy
552  double e_norm = eng / m_pivot.value();
553  double plaw = m_mc_norm * std::pow(e_norm, m_mc_exponent-1.0);
554 
555  // Compute logparabola at given energy
556  double logparabola = prefactor() *
557  std::pow(e_norm,index()+curvature()*std::log(e_norm));
558 
559  // Compute acceptance fraction
560  acceptance_fraction = logparabola / plaw;
561 
562  } while (ran.uniform() > acceptance_fraction);
563 
564  // Set energy
565  energy.MeV(eng);
566 
567  // Return energy
568  return energy;
569 }
570 
571 
572 /***********************************************************************//**
573  * @brief Read model from XML element
574  *
575  * @param[in] xml XML element.
576  *
577  * Read the spectral log parabola information from an XML element.
578  ***************************************************************************/
580 {
581  // Get parameter pointers
586 
587  // Read parameters
588  m_norm.read(*norm);
589  m_index.read(*index);
590  m_curvature.read(*curvature);
591  m_pivot.read(*pivot);
592 
593  // In case that legacy parameters were used we need to negate the index
594  // and curvature because they are differently defined
595  if (gammalib::xml_has_par(xml, "alpha") &&
596  gammalib::xml_has_par(xml, "beta")) {
597  m_index.min(-m_index.max());
598  m_index.max(-m_index.min());
603  }
604 
605  // Return
606  return;
607 }
608 
609 
610 /***********************************************************************//**
611  * @brief Write model into XML element
612  *
613  * @param[in] xml XML element.
614  *
615  * @exception GException::model_invalid_spectral
616  * Existing XML element is not of type "LogParabola"
617  *
618  * Write the LogParabola model information into an XML element.
619  ***************************************************************************/
621 {
622  // Set model type
623  if (xml.attribute("type") == "") {
624  xml.attribute("type", type());
625  }
626 
627  // Verify model type
628  if (xml.attribute("type") != type()) {
630  "Spectral model is not of type \""+type()+"\".");
631  }
632 
633  // Check if parameters are Fermi/LAT style and need to be negated
634  GModelPar inx = m_index;
635  GModelPar crv = m_curvature;
636  if (inx.name() == "alpha") {
637  inx.min(-inx.max());
638  inx.max(-inx.min());
639  inx.value(-inx.value());
640  }
641  if (crv.name() == "beta") {
642  crv.min(-crv.max());
643  crv.max(-crv.min());
644  crv.value(-crv.value());
645  }
646 
647  // Get XML parameters
652 
653  // Write parameters
654  m_norm.write(*norm);
655  inx.write(*index);
656  crv.write(*curvature);
657  m_pivot.write(*pivot);
658 
659  // Return
660  return;
661 }
662 
663 
664 /***********************************************************************//**
665  * @brief Print LogParabola information
666  *
667  * @param[in] chatter Chattiness (defaults to NORMAL).
668  * @return String containing LogParabola information.
669  ***************************************************************************/
670 std::string GModelSpectralLogParabola::print(const GChatter& chatter) const
671 {
672  // Initialise result string
673  std::string result;
674 
675  // Continue only if chatter is not silent
676  if (chatter != SILENT) {
677 
678  // Append header
679  result.append("=== GModelSpectralLogParabola ===");
680 
681  // Append information
682  result.append("\n"+gammalib::parformat("Number of parameters"));
683  result.append(gammalib::str(size()));
684  for (int i = 0; i < size(); ++i) {
685  result.append("\n"+m_pars[i]->print(chatter));
686  }
687 
688  } // endif: chatter was not silent
689 
690  // Return result
691  return result;
692 }
693 
694 
695 /*==========================================================================
696  = =
697  = Private methods =
698  = =
699  ==========================================================================*/
700 
701 /***********************************************************************//**
702  * @brief Initialise class members
703  ***************************************************************************/
705 {
706  // Initialise model type
707  m_type = "LogParabola";
708 
709  // Initialise powerlaw normalisation
710  m_norm.clear();
711  m_norm.name("Prefactor");
712  m_norm.unit("ph/cm2/s/MeV");
713  m_norm.scale(1.0);
714  m_norm.value(1.0); // default: 1.0
715  m_norm.min(0.0); // min: 0.0
716  m_norm.free();
717  m_norm.gradient(0.0);
718  m_norm.has_grad(true);
719 
720  // Initialise powerlaw index
721  m_index.clear();
722  m_index.name("Index");
723  m_index.scale(1.0);
724  m_index.value(-2.0); // default: -2.0
725  m_index.range(-10.0,+10.0); // range: [-10,+10]
726  m_index.free();
727  m_index.gradient(0.0);
728  m_index.has_grad(true);
729 
730  // Initialise curvature
731  m_curvature.clear();
732  m_curvature.name("Curvature");
733  m_curvature.scale(1.0);
734  m_curvature.value(-0.1); // default: -2.0
735  m_curvature.range(-10.0,+10.0); // range: [-10,+10]
736  m_curvature.free();
737  m_curvature.gradient(0.0);
738  m_curvature.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); // default: 100
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_curvature);
755  m_pars.push_back(&m_pivot);
756 
757  // Initialise eval cache
759  m_last_index = 1.0e30;
760  m_last_curvature = 1.0e30;
761  m_last_pivot = 1.0e30;
762  m_last_e_norm = 0.0;
763  m_last_log_e_norm = 0.0;
764  m_last_exponent = 0.0;
765  m_last_power = 0.0;
766 
767  // Initialise MC cache
768  m_mc_emin = 0.0;
769  m_mc_emax = 0.0;
770  m_mc_exponent = 0.0;
771  m_mc_pow_emin = 0.0;
772  m_mc_pow_ewidth = 0.0;
773  m_mc_norm = 0.0;
774 
775  // Return
776  return;
777 }
778 
779 
780 /***********************************************************************//**
781  * @brief Copy class members
782  *
783  * @param[in] model GModelSpectralLogParabola members which should be copied.
784  ***************************************************************************/
786 {
787  // Copy members
788  m_type = model.m_type;
789  m_norm = model.m_norm;
790  m_index = model.m_index;
791  m_curvature = model.m_curvature;
792  m_pivot = model.m_pivot;
793 
794  // Set parameter pointer(s)
795  m_pars.clear();
796  m_pars.push_back(&m_norm);
797  m_pars.push_back(&m_index);
798  m_pars.push_back(&m_curvature);
799  m_pars.push_back(&m_pivot);
800 
801  // Copy eval cache
803  m_last_index = model.m_last_index;
805  m_last_pivot = model.m_last_pivot;
809  m_last_power = model.m_last_power;
810 
811  // Copy MC cache
812  m_mc_emin = model.m_mc_emin;
813  m_mc_emax = model.m_mc_emax;
817  m_mc_norm = model.m_mc_norm;
818 
819  // Return
820  return;
821 }
822 
823 
824 /***********************************************************************//**
825  * @brief Delete class members
826  ***************************************************************************/
828 {
829  // Return
830  return;
831 }
832 
833 
834 /***********************************************************************//**
835  * @brief Update eval precomputation cache
836  *
837  * @param[in] energy Energy.
838  *
839  * Updates the precomputation cache for eval() method.
840  ***************************************************************************/
842 {
843  // Get parameter values (takes 3 multiplications which are difficult
844  // to avoid)
845  double index = m_index.value();
846  double curvature = m_curvature.value();
847  double pivot = m_pivot.value();
848 
849  // If the energy or one of the parameters index, curvature or pivot energy
850  // has changed then recompute the cache
851  if ((m_last_energy != energy) ||
852  (m_last_index != index) ||
853  (m_last_curvature != curvature) ||
854  (m_last_pivot != pivot)) {
855 
856  // Store actual energy and parameter values
857  m_last_energy = energy;
861 
862  // Compute and store value
863  double eng = energy.MeV();
864  m_last_e_norm = eng / m_last_pivot;
868 
869  } // endif: recomputation was required
870 
871  // Return
872  return;
873 }
874 
875 
876 /***********************************************************************//**
877  * @brief Update Monte Carlo pre computation cache
878  *
879  * @param[in] emin Minimum photon energy.
880  * @param[in] emax Maximum photon energy.
881  * @param[in] time True photon arrival time.
882  *
883  * Updates the precomputation cache for Monte Carlo simulations.
884  ***************************************************************************/
886  const GEnergy& emax,
887  const GTime& time) const
888 
889 {
890  // Only update if boundaries have changed
891  if(emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax){
892 
893  // Store energy boundaries
894  m_mc_emin = emin.MeV();
895  m_mc_emax = emax.MeV();
896 
897  // Find a corresponding power law with the criterion
898  // Plaw > LogParabola in the given interval
899  double index_pl;
900 
901  // Checking the sign of spectrum curvature
902  if (curvature() < 0) {
903 
904  // Use the spectral index at the pivot energy of the LogParabola
905  index_pl = index();
906  m_mc_norm = prefactor();
907  }
908  else {
909  // Use a power law which connects the ends of the convex,
910  // curved model
911 
912  // Plaw index defined by the slope of a straight line in the
913  // log-log-plane
914  index_pl = std::log(eval(emin,time)/eval(emax,time)) /
915  std::log(emin.MeV()/emax.MeV());
916 
917  // Plaw norm defined such that Plaw = LogParabola at emin
918  m_mc_norm = eval(emin,time) / std::pow(emin.MeV()/m_pivot.value(),index_pl);
919 
920  }
921 
922  // Set precomputation cache
923  if (index_pl != -1.0) {
924  m_mc_exponent = index_pl + 1.0;
927  }
928  else {
929  m_mc_exponent = 0.0;
930  m_mc_pow_emin = std::log(emin.MeV());
932  }
933 
934  }
935 
936  // Return
937  return;
938 }
const double & factor_gradient(void) const
Return parameter gradient factor.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
Abstract spectral model base class.
GModelSpectralLogParabola(void)
Void constructor.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
double m_mc_exponent
Exponent (index+1)
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
GModelPar m_index
Spectral index.
virtual void clear(void)
Clear log parabola model.
void copy_members(const GModelSpectralLogParabola &model)
Copy class members.
std::vector< GModelPar * > m_pars
Parameter pointers.
virtual GModelSpectralLogParabola * clone(void) const
Clone log parabola model.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print LogParabola information.
XML element node class.
Definition: GXmlElement.hpp:47
Spectral model registry class definition.
double max(void) const
Return parameter maximum boundary.
void update_eval_cache(const GEnergy &energy) const
Update eval precomputation cache.
GEnergy pivot(void) const
Return pivot energy.
Random number generator class.
Definition: GRan.hpp:44
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
LogParabola spectral model class.
double m_last_log_e_norm
Last ln(E/Epivot) value.
double min(void) const
Return parameter minimum boundary.
double m_last_curvature
Last curvature parameters.
#define G_WRITE
#define G_READ
bool is_free(void) const
Signal if parameter is free.
virtual std::string type(void) const
Return model type.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
double prefactor(void) const
Return pre factor.
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
Model parameter class.
Definition: GModelPar.hpp:87
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
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
double curvature(void) const
Return spectral curvature.
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
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
void clear(void)
Clear parameter.
void init_members(void)
Initialise class members.
GModelPar m_norm
Normalization factor.
virtual void read(const GXmlElement &xml)
Read model from XML element.
GChatter
Definition: GTypemaps.hpp:33
void free_members(void)
Delete class members.
double m_last_pivot
Last pivot parameter.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns Monte Carlo energy between [emin, emax].
Interface definition for the spectral model registry class.
void autoscale(void)
Autoscale parameters.
double m_mc_pow_emin
Power of minimum energy.
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1475
double index(void) const
Return spectral index.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate model.
double m_last_e_norm
Last E/Epivot value.
void eps(const double &eps)
Set relative precision.
Definition: GIntegral.hpp:206
#define G_MC
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
void init_members(void)
Initialise class members.
GEnergy m_last_energy
Last energy value.
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 m_last_index
Last index parameter.
double value(void) const
Return parameter value.
const std::string & unit(void) const
Return parameter unit.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1332
double m_mc_pow_ewidth
Power of energy width.
Exception handler interface definition.
Log parabola spectral model class definition.
virtual ~GModelSpectralLogParabola(void)
Destructor.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
const GModelSpectralLogParabola g_spectral_logp_seed1("LogParabola","Prefactor","Index","PivotEnergy","Curvature")
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void update_mc_cache(const GEnergy &emin, const GEnergy &emax, const GTime &time) const
Update Monte Carlo pre computation cache.
Integration class interface definition.
virtual GModelSpectralLogParabola & operator=(const GModelSpectralLogParabola &model)
Assignment operator.
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
double m_last_power
Last power value.
double m_mc_norm
Norm of powerlaw model at logparabola pivot energy.
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