GammaLib  1.7.0.dev
 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-2017 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())
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())
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  * @exception GException::erange_invalid
556  * Energy range is invalid (emin < emax required).
557  *
558  * Returns Monte Carlo energy by randomly drawing from a broken power law.
559  ***************************************************************************/
561  const GEnergy& emax,
562  const GTime& time,
563  GRan& ran) const
564 {
565  // Throw an exception if energy range is invalid
566  if (emin >= emax) {
567  throw GException::erange_invalid(G_MC, emin.MeV(), emax.MeV(),
568  "Minimum energy < maximum energy required.");
569  }
570 
571  // Allocate energy
572  GEnergy energy;
573 
574  // Update cache
575  update_mc_cache(emin, emax);
576 
577  // Determine in which bin we reside
578  int inx = 0;
579  if (m_mc_cum.size() > 1) {
580  double u = ran.uniform();
581  for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
582  if (m_mc_cum[inx-1] <= u)
583  break;
584  }
585  }
586 
587  // Get random energy for specific bin
588  if (m_mc_exp[inx] != 0.0) {
589  double e_min = m_mc_min[inx];
590  double e_max = m_mc_max[inx];
591  double u = ran.uniform();
592  double eng = (u > 0.0)
593  ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
594  : 0.0;
595  energy.MeV(eng);
596  }
597  else {
598  double e_min = m_mc_min[inx];
599  double e_max = m_mc_max[inx];
600  double u = ran.uniform();
601  double eng = std::exp(u * (e_max - e_min) + e_min);
602  energy.MeV(eng);
603  }
604 
605  // Return energy
606  return energy;
607 }
608 
609 
610 /***********************************************************************//**
611  * @brief Read model from XML element
612  *
613  * @param[in] xml XML element.
614  *
615  * Reads the spectral information from an XML element.
616  ***************************************************************************/
618 {
619  // Get remaining XML parameters
624 
625  // Read parameters
626  m_norm.read(*prefactor);
627  m_index1.read(*index1);
628  m_breakenergy.read(*breakenergy);
629  m_index2.read(*index2);
630 
631  // Return
632  return;
633 }
634 
635 
636 /***********************************************************************//**
637  * @brief Write model into XML element
638  *
639  * @param[in] xml XML element.
640  *
641  * @exception GException::model_invalid_spectral
642  * Existing XML element is not of the expected type.
643  *
644  * Writes the spectral information into an XML element.
645  ***************************************************************************/
647 {
648  // Set model type
649  if (xml.attribute("type") == "") {
650  xml.attribute("type", type());
651  }
652 
653  // Verify model type
654  if (xml.attribute("type") != type()) {
656  "Spectral model is not of type \""+type()+"\".");
657  }
658 
659  // Get XML parameters
664 
665  // Write parameters
666  m_norm.write(*norm);
667  m_index1.write(*index1);
668  m_index2.write(*index2);
669  m_breakenergy.write(*breakenergy);
670 
671  // Return
672  return;
673 }
674 
675 
676 /***********************************************************************//**
677  * @brief Print brokenpowerlaw information
678  *
679  * @param[in] chatter Chattiness (defaults to NORMAL).
680  * @return String containing model information.
681  ***************************************************************************/
682 std::string GModelSpectralBrokenPlaw::print(const GChatter& chatter) const
683 {
684  // Initialise result string
685  std::string result;
686 
687  // Continue only if chatter is not silent
688  if (chatter != SILENT) {
689 
690  // Append header
691  result.append("=== GModelSpectralBrokenPlaw ===");
692 
693  // Append information
694  result.append("\n"+gammalib::parformat("Number of parameters"));
695  result.append(gammalib::str(size()));
696  for (int i = 0; i < size(); ++i) {
697  result.append("\n"+m_pars[i]->print(chatter));
698  }
699 
700  } // endif: chatter was not silent
701 
702  // Return result
703  return result;
704 }
705 
706 
707 /*==========================================================================
708  = =
709  = Private methods =
710  = =
711  ==========================================================================*/
712 
713 /***********************************************************************//**
714  * @brief Initialise class members
715  ***************************************************************************/
717 {
718  // Initialise model type
719  m_type = "BrokenPowerLaw";
720 
721  // Initialise powerlaw normalisation
722  m_norm.clear();
723  m_norm.name("Prefactor");
724  m_norm.unit("ph/cm2/s/MeV");
725  m_norm.scale(1.0);
726  m_norm.value(1.0); // default: 1.0
727  m_norm.min(0.0); // min: 0.0
728  m_norm.free();
729  m_norm.gradient(0.0);
730  m_norm.has_grad(true);
731 
732  // Initialise powerlaw index1
733  m_index1.clear();
734  m_index1.name("Index1");
735  m_index1.scale(1.0);
736  m_index1.value(-2.0); // default: -2.0
737  m_index1.range(-10.0,+10.0); // range: [-10,+10]
738  m_index1.free();
739  m_index1.gradient(0.0);
740  m_index1.has_grad(true);
741 
742  // Initialise powerlaw index2
743  m_index2.clear();
744  m_index2.name("Index2");
745  m_index2.scale(1.0);
746  m_index2.value(-2.0); // default: -2.0
747  m_index2.range(-10.0,+10.0); // range: [-10,+10]
748  m_index2.free();
749  m_index2.gradient(0.0);
750  m_index2.has_grad(true);
751 
752  // Initialise break energy
754  m_breakenergy.name("BreakEnergy");
755  m_breakenergy.unit("MeV");
756  m_breakenergy.scale(1.0);
757  m_breakenergy.value(100.0); // default: 100
758  m_breakenergy.fix();
759  m_breakenergy.gradient(0.0);
760  m_breakenergy.has_grad(true);
761 
762  // Set parameter pointer(s)
763  m_pars.clear();
764  m_pars.push_back(&m_norm);
765  m_pars.push_back(&m_index1);
766  m_pars.push_back(&m_breakenergy);
767  m_pars.push_back(&m_index2);
768 
769  // Initialise eval cache
771  m_last_index1 = 1.0e30;
772  m_last_index2 = 1.0e30;
773  m_last_breakenergy = 1.0e30;
774  m_last_e_norm = 0.0;
775  m_last_log_e_norm = 0.0;
776  m_last_power = 0.0;
777 
778  // Initialise MC cache
779  m_mc_exponent1 = 0.0;
780  m_mc_exponent2 = 0.0;
781  m_mc_pow_emin = 0.0;
782  m_mc_pow_ewidth = 0.0;
783 
784  // Initialise MC cache
785  m_mc_emin = 0.0;
786  m_mc_emax = 0.0;
787 
788 
789  // Return
790  return;
791 }
792 
793 
794 /***********************************************************************//**
795  * @brief Copy class members
796  *
797  * @param[in] model GModelSpectralBrokenPlaw members which should be copied.
798  ***************************************************************************/
800 {
801  // Copy members
802  m_type = model.m_type;
803  m_norm = model.m_norm;
804  m_index1 = model.m_index1;
806  m_index2 = model.m_index2;
807 
808  // Set parameter pointer(s)
809  m_pars.clear();
810  m_pars.push_back(&m_norm);
811  m_pars.push_back(&m_index1);
812  m_pars.push_back(&m_breakenergy);
813  m_pars.push_back(&m_index2);
814 
815  // Copy eval cache
822  m_last_power = model.m_last_power;
823 
824  // Copy MC cache
825  m_mc_emin = model.m_mc_emin;
826  m_mc_emax = model.m_mc_emax;
831 
832  // Return
833  return;
834 }
835 
836 
837 /***********************************************************************//**
838  * @brief Delete class members
839  ***************************************************************************/
841 {
842  // Return
843  return;
844 }
845 
846 
847 /***********************************************************************//**
848  * @brief Update eval precomputation cache
849  *
850  * @param[in] energy Energy.
851  *
852  * Updates the precomputation cache for eval() the method.
853  ***************************************************************************/
855 {
856  // Get parameter values (takes 2 multiplications which are difficult
857  // to avoid)
858  double index1 = m_index1.value();
859  double breakenergy = m_breakenergy.value();
860  double index2 = m_index2.value();
861 
862  // If the energy or one of the parameters index1, index2 or breakenergy
863  // energy has changed then recompute the cache
864  if ((m_last_energy != energy) ||
865  (m_last_index1 != index1) ||
866  (m_last_index2 != index2) ||
867  (m_last_breakenergy != breakenergy)) {
868 
869  // Store actual energy and parameter values
870  m_last_energy = energy;
874 
875  // Compute and store value
876  double eng = energy.MeV();
879  if (eng < m_last_breakenergy ) {
881  }
882  else {
884  }
885 
886  } // endif: recomputation was required
887 
888  // Return
889  return;
890 }
891 
892 
893 /***********************************************************************//**
894  * @brief Update Monte Carlo pre computation cache
895  *
896  * @param[in] emin Minimum photon energy.
897  * @param[in] emax Maximum photon energy.
898  *
899  * Updates the precomputation cache for Monte Carlo simulations.
900  ***************************************************************************/
902  const GEnergy& emax) const
903 {
904  // Check if we need to update the cache
905  if (emin.MeV() != m_mc_emin || emax.MeV() != m_mc_emax) {
906 
907  // Store new energy interval
908  m_mc_emin = emin.MeV();
909  m_mc_emax = emax.MeV();
910 
911  // Initialise cache
912  m_mc_cum.clear();
913  m_mc_min.clear();
914  m_mc_max.clear();
915  m_mc_exp.clear();
916 
917  // Get energy range in MeV
918  double e_min = emin.MeV();
919  double e_max = emax.MeV();
920 
921  // Continue only if e_max > e_min
922  if (e_max > e_min) {
923 
924  // Allocate flux
925  double flux;
926 
927  // Determine left node index for minimum and maximum energy
928  int inx_emin = (e_min < m_breakenergy.value() ) ? 0 : 1;
929  int inx_emax = (e_max < m_breakenergy.value() ) ? 0 : 1;
930 
931  // If both energies are within the same node then just
932  // add this one node on the stack
933  if (inx_emin == inx_emax) {
934  double exp_valid = (e_min < m_breakenergy.value())
935  ? m_index1.value() : m_index2.value();
936  flux = m_norm.value() *
938  e_max,
940  exp_valid);
941  m_mc_cum.push_back(flux);
942  m_mc_min.push_back(e_min);
943  m_mc_max.push_back(e_max);
944  m_mc_exp.push_back(exp_valid);
945  }
946 
947  // ... otherwise integrate over both nodes
948  else {
949  // just enter the values for first pl: bin [0]
950  flux = m_norm.value() *
954  m_index1.value());
955  m_mc_cum.push_back(flux);
956  m_mc_exp.push_back(m_index1.value());
957  m_mc_min.push_back(e_min);
958  m_mc_max.push_back(m_breakenergy.value());
959 
960  // and for
961  flux = m_norm.value() *
963  e_max,
965  m_index2.value());
966  m_mc_cum.push_back(flux);
967  m_mc_exp.push_back(m_index2.value());
968  m_mc_max.push_back(e_max);
969  m_mc_min.push_back(m_breakenergy.value());
970  } // endelse: emin and emax not between same nodes
971 
972  // Build cumulative distribution
973  for (int i = 1; i < m_mc_cum.size(); ++i) {
974  m_mc_cum[i] += m_mc_cum[i-1];
975  }
976  double norm = m_mc_cum[m_mc_cum.size()-1];
977  for (int i = 0; i < m_mc_cum.size(); ++i) {
978  m_mc_cum[i] /= norm;
979  }
980 
981  // Set MC values
982  for (int i = 0; i < m_mc_cum.size(); ++i) {
983 
984  // Compute exponent
985  double exponent = m_mc_exp[i] + 1.0;
986 
987  // Exponent dependend computation
988  if (std::abs(exponent) > 1.0e-11) {
989  m_mc_exp[i] = exponent;
990  m_mc_min[i] = std::pow(m_mc_min[i], exponent);
991  m_mc_max[i] = std::pow(m_mc_max[i], exponent);
992  }
993  else {
994  m_mc_exp[i] = 0.0;
995  m_mc_min[i] = std::log(m_mc_min[i]);
996  m_mc_max[i] = std::log(m_mc_max[i]);
997  }
998 
999  } // endfor: set MC values
1000 
1001  } // endif: e_max > e_min
1002 
1003  } // endif: Update was required
1004 
1005  // Return
1006  return;
1007 }
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 gradient factor.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:821
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:1163
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:353
virtual std::string type(void) const
Return model type.
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:47
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:54
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:1127
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:184
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:167
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.
const GXmlAttribute * attribute(const int &index) const
Return attribute.
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:1513
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:1184
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
void clear(void)
Clear parameter.
double 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
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:1332
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:1142
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:1022
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:1562
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 value factor.
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:1079
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413