GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralBrokenPlaw.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralBrokenPlaw.cpp - Broken power law spectrum class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2013-2021 by Anneli Schulz *
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 GModelSpectralBrokenPlaw.cpp
23  * @brief Broken power law spectrum class implementation
24  * @author Anneli Schulz
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"
37 
38 /* __ Constants __________________________________________________________ */
39 
40 /* __ Globals ____________________________________________________________ */
41 const GModelSpectralBrokenPlaw g_spectral_blaw_seed1("BrokenPowerLaw",
42  "Prefactor",
43  "Index1",
44  "BreakEnergy",
45  "Index2");
46 const GModelSpectralRegistry g_spectral_blaw_registry1(&g_spectral_blaw_seed1);
47 #if defined(G_LEGACY_XML_FORMAT)
48 const GModelSpectralBrokenPlaw g_spectral_blaw_seed2("BrokenPowerLaw",
49  "Prefactor",
50  "Index1",
51  "BreakValue",
52  "Index2");
53 const GModelSpectralRegistry g_spectral_blaw_registry2(&g_spectral_blaw_seed2);
54 #endif
55 
56 /* __ Method name definitions ____________________________________________ */
57 #define G_MC "GModelSpectralBrokenPlaw::mc(GEnergy&, GEnergy&, GTime&,"\
58  " GRan&)"
59 #define G_READ "GModelSpectralBrokenPlaw::read(GXmlElement&)"
60 #define G_WRITE "GModelSpectralBrokenPlaw::write(GXmlElement&)"
61 
62 /* __ Macros _____________________________________________________________ */
63 
64 /* __ Coding definitions _________________________________________________ */
65 
66 /* __ Debug definitions __________________________________________________ */
67 
68 
69 /*==========================================================================
70  = =
71  = Constructors/destructors =
72  = =
73  ==========================================================================*/
74 
75 /***********************************************************************//**
76  * @brief Void constructor
77  ***************************************************************************/
79 {
80  // Initialise private members for clean destruction
81  init_members();
82 
83  // Return
84  return;
85 }
86 
87 
88 /***********************************************************************//**
89  * @brief Model type and parameter name constructor
90  *
91  * @param[in] type Model type.
92  * @param[in] prefactor Name of prefactor parameter.
93  * @param[in] index1 Name of index1 parameter.
94  * @param[in] breakenergy Name of breakenergy parameter.
95  * @param[in] index2 Name of index2 parameter.
96  ***************************************************************************/
98  const std::string& prefactor,
99  const std::string& index1,
100  const std::string& breakenergy,
101  const std::string& index2) :
103 {
104  // Initialise members
105  init_members();
106 
107  // Set model type
108  m_type = type;
109 
110  // Set parameter names
111  m_norm.name(prefactor);
112  m_index1.name(index1);
113  m_breakenergy.name(breakenergy);
114  m_index2.name(index2);
115 
116  // Return
117  return;
118 }
119 
120 
121 /***********************************************************************//**
122  * @brief Parameter constructor
123  *
124  * @param[in] prefactor Power law pre factor (ph/cm2/s/MeV).
125  * @param[in] index1 Power law index1.
126  * @param[in] breakenergy break energy.
127  * @param[in] index2 Power law index1.
128  *
129  * Constructs a spectral power law using the model parameters
130  *
131  * power law @p prefactor (ph/cm2/s/MeV)
132  * spectral @p index1
133  * Energy @p breakenergy of spectral break
134  * spectral @p index2
135  ***************************************************************************/
137  const double& index1,
138  const GEnergy& breakenergy,
139  const double& index2) :
141 {
142  // Initialise members
143  init_members();
144 
145  // Set parameters
146  m_norm.value(prefactor);
147  m_index1.value(index1);
148  m_breakenergy.value(breakenergy.MeV()); // Internally stored in MeV
149  m_index2.value(index2);
150 
151  // Perform autoscaling of parameter
152  autoscale();
153 
154  // Return
155  return;
156 }
157 
158 
159 /***********************************************************************//**
160  * @brief XML constructor
161  *
162  * @param[in] xml XML element containing position information.
163  *
164  * Constructs a power law spectral model by extracting information from an
165  * XML element. See the read() method for more information about the expected
166  * structure of the XML element.
167  ***************************************************************************/
170 {
171  // Initialise members
172  init_members();
173 
174  // Read information from XML element
175  read(xml);
176 
177  // Return
178  return;
179 }
180 
181 
182 /***********************************************************************//**
183  * @brief Copy constructor
184  *
185  * @param[in] model Broken power law model.
186  ***************************************************************************/
188  GModelSpectral(model)
189 {
190  // Initialise members
191  init_members();
192 
193  // Copy members
194  copy_members(model);
195 
196  // Return
197  return;
198 }
199 
200 
201 /***********************************************************************//**
202  * @brief Destructor
203  ***************************************************************************/
205 {
206  // Free members
207  free_members();
208 
209  // Return
210  return;
211 }
212 
213 
214 /*==========================================================================
215  = =
216  = Operators =
217  = =
218  ==========================================================================*/
219 
220 /***********************************************************************//**
221  * @brief Assignment operator
222  *
223  * @param[in] model Broken power law model.
224  * @return Broken power law model.
225  ***************************************************************************/
227 {
228  // Execute only if object is not identical
229  if (this != &model) {
230 
231  // Copy base class members
232  this->GModelSpectral::operator=(model);
233 
234  // Free members
235  free_members();
236 
237  // Initialise members
238  init_members();
239 
240  // Copy members
241  copy_members(model);
242 
243  } // endif: object was not identical
244 
245  // Return
246  return *this;
247 }
248 
249 
250 /*==========================================================================
251  = =
252  = Public methods =
253  = =
254  ==========================================================================*/
255 
256 /***********************************************************************//**
257  * @brief Clear broken power law model
258  ***************************************************************************/
260 {
261  // Free class members (base and derived classes, derived class first)
262  free_members();
264 
265  // Initialise members
267  init_members();
268 
269  // Return
270  return;
271 }
272 
273 
274 /***********************************************************************//**
275  * @brief Clone broken power law model
276  *
277  * @return Pointer to deep copy of broken power law model.
278  ***************************************************************************/
280 {
281  // Clone spectral broken power law model
282  return new GModelSpectralBrokenPlaw(*this);
283 }
284 
285 
286 /***********************************************************************//**
287  * @brief Evaluate function
288  *
289  * @param[in] srcEng True photon energy.
290  * @param[in] srcTime True photon arrival time.
291  * @param[in] gradients Compute gradients?
292  * @return Model value (ph/cm2/s/MeV).
293  *
294  * Evaluates
295  *
296  * \f[
297  * S_{\rm E}(E | t) = k_0 \times \left \{
298  * \begin{eqnarray}
299  * \left( \frac{E}{E_b} \right)^{\gamma_1} & {\rm if\,\,} E < E_b \\
300  * \left( \frac{E}{E_b} \right)^{\gamma_2} & {\rm otherwise}
301  * \end{eqnarray}
302  * \right .
303  * \f]
304  *
305  * where
306  * \f$k_0\f$ is the normalization or prefactor,
307  * \f$\gamma_1\f$ is the spectral index before the break,
308  * \f$\gamma_2\f$ is the spectral index after the break, and
309  * \f$E_b\f$ is the break energy.
310  *
311  * If the @p gradients flag is true the method will also compute the
312  * partial derivatives of the model with respect to the parameters using
313  *
314  * \f[
315  * \frac{\delta S_{\rm E}(E | t)}{\delta k_0} =
316  * \frac{S_{\rm E}(E | t)}{k_0}
317  * \f]
318  *
319  * \f[
320  * \frac{\delta S_{\rm E}(E | t)}{\delta \gamma_1} = \left \{
321  * \begin{eqnarray}
322  * S_{\rm E}(E | t) \, \ln(E/E_b) & {\rm if\,\,} E < E_b \\
323  * 0 & {\rm otherwise}
324  * \end{eqnarray}
325  * \right .
326  * \f]
327  *
328  * \f[
329  * \frac{\delta S_{\rm E}(E | t)}{\delta \gamma_2} = \left \{
330  * \begin{eqnarray}
331  * 0 & {\rm if\,\,} E < E_b \\
332  * S_{\rm E}(E | t) \, \ln(E/E_b) & {\rm otherwise}
333  * \end{eqnarray}
334  * \right .
335  * \f]
336  *
337  * \f[
338  * \frac{\delta S_{\rm E}(E | t)}{\delta E_b} = k_0 \times \left \{
339  * \begin{eqnarray}
340  * -S_{\rm E}(E | t) \, \left( \frac{\gamma_1}{E_b} \right) & {\rm if\,\,} E < E_b \\
341  * -S_{\rm E}(E | t) \, \left( \frac{\gamma_2}{E_b} \right) & {\rm otherwise}
342  * \end{eqnarray}
343  * \right .
344  * \f]
345  *
346  * @todo The method expects that energy!=0. Otherwise Inf or NaN may result.
347  ***************************************************************************/
349  const GTime& srcTime,
350  const bool& gradients) const
351 {
352  // Update the evaluation cache
353  update_eval_cache(srcEng);
354 
355  // Compute function value
356  double value = m_norm.value() * m_last_power;
357 
358  // Optionally compute gradients
359  if (gradients) {
360 
361  // Compute normalisation gradient
362  double g_norm = (m_norm.is_free()) ? m_norm.scale() * m_last_power : 0.0;
363  m_norm.factor_gradient(g_norm);
364 
365  // Compute index and break value gradients
366  if (srcEng.MeV() < m_breakenergy.value()) {
367  double g_index = (m_index1.is_free())
368  ? value * m_index1.scale() * m_last_log_e_norm
369  : 0.0;
370  double g_break = (m_breakenergy.is_free() && m_breakenergy.factor_value() != 0.0)
372  : 0.0;
373  m_index1.factor_gradient(g_index);
376  }
377  else {
378  double g_index = (m_index2.is_free())
379  ? value * m_index2.scale() * m_last_log_e_norm
380  : 0.0;
381  double g_break = (m_breakenergy.is_free() && m_breakenergy.factor_value() != 0.0)
383  : 0.0;
385  m_index2.factor_gradient(g_index);
387  }
388 
389  } // endif: gradient computation was requested
390 
391  // Compile option: Check for NaN/Inf
392  #if defined(G_NAN_CHECK)
393  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
394  std::cout << "*** ERROR: GModelSpectralBrokenPlaw::eval";
395  std::cout << "(srcEng=" << srcEng;
396  std::cout << ", srcTime=" << srcTime << "):";
397  std::cout << " NaN/Inf encountered";
398  std::cout << " (value=" << value;
399  std::cout << ", m_norm=" << m_norm.value();
400  std::cout << ", m_index1=" << m_index1.value();
401  std::cout << ", m_breakenergy=" << m_breakenergy.value();
402  std::cout << ", m_index2=" << m_index2.value();
403  std::cout << ", m_last_power=" << m_last_power;
404  std::cout << ")" << std::endl;
405  }
406  #endif
407 
408  // Return
409  return value;
410 }
411 
412 
413 /***********************************************************************//**
414  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
415  *
416  * @param[in] emin Minimum photon energy.
417  * @param[in] emax Maximum photon energy.
418  * @return Photon flux (ph/cm2/s).
419  *
420  * Computes
421  *
422  * \f[
423  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
424  * \f]
425  *
426  * where
427  * - [@p emin, @p emax] is an energy interval, and
428  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
429  * The integration is done analytically.
430  ***************************************************************************/
432  const GEnergy& emax) const
433 {
434  // Initialise flux
435  double flux = 0.0;
436 
437  // Compute only if integration range is valid
438  if (emin < emax) {
439 
440  // First case: emin < breakenergy < emax
441  if (emin.MeV() < m_breakenergy.value() && m_breakenergy.value() < emax.MeV()) {
442 
443  // Compute photon flux
444  flux = m_norm.value() *
448  m_index1.value()) +
450  emax.MeV(),
452  m_index2.value()));
453  }
454 
455  // Second case: breakenergy >= emax
456  else if (m_breakenergy.value() >= emax.MeV()) {
457  flux = m_norm.value() *
459  emax.MeV(),
461  m_index1.value());
462  }
463 
464  // Third case: breakenergy <= emin
465  else if (m_breakenergy.value() <= emin.MeV()) {
466  flux = m_norm.value() *
468  emax.MeV(),
470  m_index2.value());
471  }
472 
473  } // endif: integration range was valid
474 
475  // Return flux
476  return flux;
477 }
478 
479 
480 /***********************************************************************//**
481  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
482  *
483  * @param[in] emin Minimum photon energy.
484  * @param[in] emax Maximum photon energy.
485  * @return Energy flux (erg/cm2/s).
486  *
487  * Computes
488  *
489  * \f[
490  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
491  * \f]
492  *
493  * where
494  * - [@p emin, @p emax] is an energy interval, and
495  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
496  * The integration is done analytically.
497  ***************************************************************************/
499  const GEnergy& emax) const
500 {
501  // Initialise flux
502  double eflux = 0.0;
503 
504  // Compute only if integration range is valid
505  if (emin < emax) {
506 
507  // First case: emin < breakenergy < emax
508  if (emin.MeV() < m_breakenergy.value() && m_breakenergy.value() < emax.MeV()) {
509  eflux = m_norm.value() *
513  m_index1.value()) +
515  emax.MeV(),
517  m_index2.value()));
518  }
519 
520  // Second case: breakenergy >= emax
521  else if (m_breakenergy.value() >= emax.MeV()) {
522  eflux = m_norm.value() *
524  emax.MeV(),
526  m_index1.value());
527  }
528 
529  // Third case: breakenergy <= emin
530  else if (m_breakenergy.value() <= emin.MeV()) {
531  eflux = m_norm.value() *
533  emax.MeV(),
535  m_index2.value());
536  }
537  eflux *= gammalib::MeV2erg;
538 
539  } // endif: integration range was valid
540 
541  // Return flux
542  return eflux;
543 }
544 
545 
546 /***********************************************************************//**
547  * @brief Returns Monte Carlo energy between [emin, emax]
548  *
549  * @param[in] emin Minimum photon energy.
550  * @param[in] emax Maximum photon energy.
551  * @param[in] time True photon arrival time.
552  * @param[in,out] ran Random number generator.
553  * @return Energy.
554  *
555  * Returns Monte Carlo energy by randomly drawing from a broken power law.
556  ***************************************************************************/
558  const GEnergy& emax,
559  const GTime& time,
560  GRan& ran) const
561 {
562  // Check energy interval
564 
565  // Allocate energy
566  GEnergy energy;
567 
568  // Update cache
569  update_mc_cache(emin, emax);
570 
571  // Determine in which bin we reside
572  int inx = 0;
573  if (m_mc_cum.size() > 1) {
574  double u = ran.uniform();
575  for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
576  if (m_mc_cum[inx-1] <= u)
577  break;
578  }
579  }
580 
581  // Get random energy for specific bin
582  if (m_mc_exp[inx] != 0.0) {
583  double e_min = m_mc_min[inx];
584  double e_max = m_mc_max[inx];
585  double u = ran.uniform();
586  double eng = (u > 0.0)
587  ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
588  : 0.0;
589  energy.MeV(eng);
590  }
591  else {
592  double e_min = m_mc_min[inx];
593  double e_max = m_mc_max[inx];
594  double u = ran.uniform();
595  double eng = std::exp(u * (e_max - e_min) + e_min);
596  energy.MeV(eng);
597  }
598 
599  // Return energy
600  return energy;
601 }
602 
603 
604 /***********************************************************************//**
605  * @brief Read model from XML element
606  *
607  * @param[in] xml XML element.
608  *
609  * Reads the spectral information from an XML element.
610  ***************************************************************************/
612 {
613  // Verify number of model parameters
615 
616  // Get remaining XML parameters
621 
622  // Read parameters
623  m_norm.read(*prefactor);
624  m_index1.read(*index1);
625  m_breakenergy.read(*breakenergy);
626  m_index2.read(*index2);
627 
628  // Return
629  return;
630 }
631 
632 
633 /***********************************************************************//**
634  * @brief Write model into XML element
635  *
636  * @param[in] xml XML element.
637  *
638  * Writes the spectral information into an XML element.
639  ***************************************************************************/
641 {
642  // Verify model type
644 
645  // Get XML parameters
650 
651  // Write parameters
652  m_norm.write(*norm);
653  m_index1.write(*index1);
654  m_index2.write(*index2);
655  m_breakenergy.write(*breakenergy);
656 
657  // Return
658  return;
659 }
660 
661 
662 /***********************************************************************//**
663  * @brief Print brokenpowerlaw information
664  *
665  * @param[in] chatter Chattiness (defaults to NORMAL).
666  * @return String containing model information.
667  ***************************************************************************/
668 std::string GModelSpectralBrokenPlaw::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("=== GModelSpectralBrokenPlaw ===");
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 = "BrokenPowerLaw";
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 index1
719  m_index1.clear();
720  m_index1.name("Index1");
721  m_index1.scale(1.0);
722  m_index1.value(-2.0); // default: -2.0
723  m_index1.range(-10.0,+10.0); // range: [-10,+10]
724  m_index1.free();
725  m_index1.gradient(0.0);
726  m_index1.has_grad(true);
727 
728  // Initialise powerlaw index2
729  m_index2.clear();
730  m_index2.name("Index2");
731  m_index2.scale(1.0);
732  m_index2.value(-2.0); // default: -2.0
733  m_index2.range(-10.0,+10.0); // range: [-10,+10]
734  m_index2.free();
735  m_index2.gradient(0.0);
736  m_index2.has_grad(true);
737 
738  // Initialise break energy
740  m_breakenergy.name("BreakEnergy");
741  m_breakenergy.unit("MeV");
742  m_breakenergy.scale(1.0);
743  m_breakenergy.value(100.0); // default: 100
744  m_breakenergy.fix();
745  m_breakenergy.gradient(0.0);
746  m_breakenergy.has_grad(true);
747 
748  // Set parameter pointer(s)
749  m_pars.clear();
750  m_pars.push_back(&m_norm);
751  m_pars.push_back(&m_index1);
752  m_pars.push_back(&m_breakenergy);
753  m_pars.push_back(&m_index2);
754 
755  // Initialise eval cache
757  m_last_index1 = 1.0e30;
758  m_last_index2 = 1.0e30;
759  m_last_breakenergy = 1.0e30;
760  m_last_e_norm = 0.0;
761  m_last_log_e_norm = 0.0;
762  m_last_power = 0.0;
763 
764  // Initialise MC cache
765  m_mc_exponent1 = 0.0;
766  m_mc_exponent2 = 0.0;
767  m_mc_pow_emin = 0.0;
768  m_mc_pow_ewidth = 0.0;
769 
770  // Initialise MC cache
771  m_mc_emin = 0.0;
772  m_mc_emax = 0.0;
773 
774 
775  // Return
776  return;
777 }
778 
779 
780 /***********************************************************************//**
781  * @brief Copy class members
782  *
783  * @param[in] model GModelSpectralBrokenPlaw members which should be copied.
784  ***************************************************************************/
786 {
787  // Copy members
788  m_type = model.m_type;
789  m_norm = model.m_norm;
790  m_index1 = model.m_index1;
792  m_index2 = model.m_index2;
793 
794  // Set parameter pointer(s)
795  m_pars.clear();
796  m_pars.push_back(&m_norm);
797  m_pars.push_back(&m_index1);
798  m_pars.push_back(&m_breakenergy);
799  m_pars.push_back(&m_index2);
800 
801  // Copy eval cache
808  m_last_power = model.m_last_power;
809 
810  // Copy MC cache
811  m_mc_emin = model.m_mc_emin;
812  m_mc_emax = model.m_mc_emax;
817 
818  // Return
819  return;
820 }
821 
822 
823 /***********************************************************************//**
824  * @brief Delete class members
825  ***************************************************************************/
827 {
828  // Return
829  return;
830 }
831 
832 
833 /***********************************************************************//**
834  * @brief Update eval precomputation cache
835  *
836  * @param[in] energy Energy.
837  *
838  * Updates the precomputation cache for eval() the method.
839  ***************************************************************************/
841 {
842  // Get parameter values (takes 2 multiplications which are difficult
843  // to avoid)
844  double index1 = m_index1.value();
845  double breakenergy = m_breakenergy.value();
846  double index2 = m_index2.value();
847 
848  // If the energy or one of the parameters index1, index2 or breakenergy
849  // energy has changed then recompute the cache
850  if ((m_last_energy != energy) ||
851  (m_last_index1 != index1) ||
852  (m_last_index2 != index2) ||
853  (m_last_breakenergy != breakenergy)) {
854 
855  // Store actual energy and parameter values
856  m_last_energy = energy;
860 
861  // Compute and store value
862  double eng = energy.MeV();
865  if (eng < m_last_breakenergy ) {
867  }
868  else {
870  }
871 
872  } // endif: recomputation was required
873 
874  // Return
875  return;
876 }
877 
878 
879 /***********************************************************************//**
880  * @brief Update Monte Carlo pre computation cache
881  *
882  * @param[in] emin Minimum photon energy.
883  * @param[in] emax Maximum photon energy.
884  *
885  * Updates the precomputation cache for Monte Carlo simulations.
886  ***************************************************************************/
888  const GEnergy& emax) const
889 {
890  // Check if we need to update the cache
891  if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) {
892 
893  // Store new energy interval
894  m_mc_emin = emin.MeV();
895  m_mc_emax = emax.MeV();
896 
897  // Initialise cache
898  m_mc_cum.clear();
899  m_mc_min.clear();
900  m_mc_max.clear();
901  m_mc_exp.clear();
902 
903  // Get energy range in MeV
904  double e_min = emin.MeV();
905  double e_max = emax.MeV();
906 
907  // Continue only if e_max > e_min
908  if (e_max > e_min) {
909 
910  // Allocate flux
911  double flux;
912 
913  // Determine left node index for minimum and maximum energy
914  int inx_emin = (e_min < m_breakenergy.value() ) ? 0 : 1;
915  int inx_emax = (e_max < m_breakenergy.value() ) ? 0 : 1;
916 
917  // If both energies are within the same node then just
918  // add this one node on the stack
919  if (inx_emin == inx_emax) {
920  double exp_valid = (e_min < m_breakenergy.value())
921  ? m_index1.value() : m_index2.value();
922  flux = m_norm.value() *
924  e_max,
926  exp_valid);
927  m_mc_cum.push_back(flux);
928  m_mc_min.push_back(e_min);
929  m_mc_max.push_back(e_max);
930  m_mc_exp.push_back(exp_valid);
931  }
932 
933  // ... otherwise integrate over both nodes
934  else {
935  // just enter the values for first pl: bin [0]
936  flux = m_norm.value() *
940  m_index1.value());
941  m_mc_cum.push_back(flux);
942  m_mc_exp.push_back(m_index1.value());
943  m_mc_min.push_back(e_min);
944  m_mc_max.push_back(m_breakenergy.value());
945 
946  // and for
947  flux = m_norm.value() *
949  e_max,
951  m_index2.value());
952  m_mc_cum.push_back(flux);
953  m_mc_exp.push_back(m_index2.value());
954  m_mc_max.push_back(e_max);
955  m_mc_min.push_back(m_breakenergy.value());
956  } // endelse: emin and emax not between same nodes
957 
958  // Build cumulative distribution
959  for (int i = 1; i < m_mc_cum.size(); ++i) {
960  m_mc_cum[i] += m_mc_cum[i-1];
961  }
962  double norm = m_mc_cum[m_mc_cum.size()-1];
963  for (int i = 0; i < m_mc_cum.size(); ++i) {
964  m_mc_cum[i] /= norm;
965  }
966 
967  // Set MC values
968  for (int i = 0; i < m_mc_cum.size(); ++i) {
969 
970  // Compute exponent
971  double exponent = m_mc_exp[i] + 1.0;
972 
973  // Exponent dependend computation
974  if (std::abs(exponent) > 1.0e-11) {
975  m_mc_exp[i] = exponent;
976  m_mc_min[i] = std::pow(m_mc_min[i], exponent);
977  m_mc_max[i] = std::pow(m_mc_max[i], exponent);
978  }
979  else {
980  m_mc_exp[i] = 0.0;
981  m_mc_min[i] = std::log(m_mc_min[i]);
982  m_mc_max[i] = std::log(m_mc_max[i]);
983  }
984 
985  } // endfor: set MC values
986 
987  } // endif: e_max > e_min
988 
989  } // endif: Update was required
990 
991  // Return
992  return;
993 }
GModelPar m_index2
Spectral index2.
void copy_members(const GModelSpectralBrokenPlaw &model)
Copy class members.
double m_last_index2
Last index1 parameter.
std::vector< double > m_mc_min
Lower boundary for MC.
const double & factor_gradient(void) const
Return parameter factor gradient.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
const std::string & name(void) const
Return parameter name.
double index2(void) const
Return power law index2.
Random number generator class definition.
void update_mc_cache(const GEnergy &emin, const GEnergy &emax) const
Update Monte Carlo pre computation cache.
Abstract spectral model base class.
virtual void read(const GXmlElement &xml)
Read model from XML element.
double gradient(void) const
Return parameter gradient.
int size(void) const
Return number of parameters.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
virtual std::string type(void) const
Return model type.
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:48
Spectral model registry class definition.
void free_members(void)
Delete class members.
#define G_MC
Random number generator class.
Definition: GRan.hpp:44
double m_last_log_e_norm
Last ln(E/Ebreakenergy) value.
#define G_WRITE
std::vector< double > m_mc_exp
Exponent for MC.
double prefactor(void) const
Return pre factor.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
GModelPar m_index1
Spectral index1.
Time class.
Definition: GTime.hpp:55
double m_last_power
Last power value.
Gammalib tools definition.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
double m_last_index1
Last index1 parameter.
virtual GModelSpectralBrokenPlaw * clone(void) const
Clone broken power law model.
void update_eval_cache(const GEnergy &energy) const
Update eval precomputation cache.
virtual void clear(void)
Clear broken power law model.
double min(void) const
Return parameter minimum boundary.
double m_mc_exponent1
Exponent (index1+1)
double plaw_energy_flux(const double &emin, const double &emax, const double &epivot, const double &gamma)
Compute energy flux between two energies for a power law.
Definition: GTools.cpp:1248
bool is_free(void) const
Signal if parameter is free.
std::vector< double > m_mc_cum
Cumulative distribution.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
const GModelSpectralBrokenPlaw g_spectral_blaw_seed1("BrokenPowerLaw","Prefactor","Index1","BreakEnergy","Index2")
GModelPar m_norm
Normalization factor.
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
virtual void write(GXmlElement &xml) const
Write model into XML element.
GEnergy breakenergy(void) const
Return breakenergy energy.
GModelSpectralBrokenPlaw(void)
Void constructor.
std::vector< double > m_mc_max
Upper boundary for MC.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
void free(void)
Free a parameter.
void fix(void)
Fix a parameter.
GModelPar m_breakenergy
Energy of spectral break.
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:1637
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns Monte Carlo energy between [emin, emax].
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void check_energy_interval(const std::string &origin, const GEnergy &emin, const GEnergy &emax)
Checks energy interval.
Definition: GException.cpp:364
void clear(void)
Clear parameter.
double m_mc_pow_emin
Power of minimum energy.
Broken power law spectrum class definition.
GChatter
Definition: GTypemaps.hpp:33
double index1(void) const
Return power law index1.
Broken power law spectral model class.
#define G_READ
void xml_check_parnum(const std::string &origin, const GXmlElement &xml, const int &number)
Checks number of parameters.
Definition: GTools.cpp:1777
Interface definition for the spectral model registry class.
void autoscale(void)
Autoscale parameters.
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
double m_last_breakenergy
Last breakenergy parameter.
virtual ~GModelSpectralBrokenPlaw(void)
Destructor.
virtual GModelSpectralBrokenPlaw & operator=(const GModelSpectralBrokenPlaw &model)
Assignment operator.
void init_members(void)
Initialise class members.
double m_mc_pow_ewidth
Power of energy width.
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.
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
const std::string & unit(void) const
Return parameter unit.
double m_mc_exponent2
Exponent (index2+1)
Exception handler interface definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
GEnergy m_last_energy
Last energy value.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
double m_last_e_norm
Last E/Ebreakenergy value.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print brokenpowerlaw information.
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:1689
void free_members(void)
Delete class members.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
void init_members(void)
Initialise class members.
const double & factor_value(void) const
Return parameter factor value.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
double plaw_photon_flux(const double &emin, const double &emax, const double &epivot, const double &gamma)
Compute photon flux between two energies for a power law.
Definition: GTools.cpp:1200
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819