GammaLib  2.1.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GModelSpectralNodes.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GModelSpectralNodes.cpp - Spectral nodes model class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2021 by Juergen Knoedlseder *
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 GModelSpectralNodes.cpp
23  * @brief Spectral nodes model class implementation
24  * @author Juergen Knoedlseder
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 "GEnergies.hpp"
36 #include "GModelSpectralNodes.hpp"
38 
39 /* __ Constants __________________________________________________________ */
40 
41 /* __ Globals ____________________________________________________________ */
43 const GModelSpectralRegistry g_spectral_nodes_registry(&g_spectral_nodes_seed);
44 
45 /* __ Method name definitions ____________________________________________ */
46 #define G_FLUX "GModelSpectralNodes::flux(GEnergy&, GEnergy&)"
47 #define G_EFLUX "GModelSpectralNodes::eflux(GEnergy&, GEnergy&)"
48 #define G_MC "GModelSpectralNodes::mc(GEnergy&, GEnergy&, GTime&, GRan&)"
49 #define G_READ "GModelSpectralNodes::read(GXmlElement&)"
50 #define G_WRITE "GModelSpectralNodes::write(GXmlElement&)"
51 #define G_APPEND "GModelSpectralNodes::append(GEnergy&, double&)"
52 #define G_INSERT "GModelSpectralNodes::insert(int&, GEnergy&, double&)"
53 #define G_REMOVE "GModelSpectralNodes::remove(int&)"
54 #define G_ENERGY_GET "GModelSpectralNodes::energy(int&)"
55 #define G_ENERGY_SET "GModelSpectralNodes::energy(int&, GEnergy&)"
56 #define G_INTENSITY_GET "GModelSpectralNodes::intensity(int&)"
57 #define G_INTENSITY_SET "GModelSpectralNodes::intensity(int&, double&)"
58 #define G_ERROR_GET "GModelSpectralNodes::error(int&)"
59 
60 /* __ Macros _____________________________________________________________ */
61 
62 /* __ Coding definitions _________________________________________________ */
63 
64 /* __ Debug definitions __________________________________________________ */
65 
66 
67 /*==========================================================================
68  = =
69  = Constructors/destructors =
70  = =
71  ==========================================================================*/
72 
73 /***********************************************************************//**
74  * @brief Void constructor
75  *
76  * Constructs an empty spectral node model.
77  ***************************************************************************/
79 {
80  // Initialise members
81  init_members();
82 
83  // Return
84  return;
85 }
86 
87 
88 /***********************************************************************//**
89  * @brief Spectral model constructor
90  *
91  * @param[in] model Spectral model.
92  * @param[in] energies Node energies.
93  *
94  * Constructs a spectral node model from any spectral model.
95  ***************************************************************************/
97  const GEnergies& energies) :
99 {
100  // Initialise members
101  init_members();
102 
103  // Reserve space for nodes
104  reserve(energies.size());
105 
106  // Append nodes for all energies
107  for (int i = 0; i < energies.size(); ++i) {
108  append(energies[i], model.eval(energies[i]));
109  }
110 
111  // Return
112  return;
113 }
114 
115 
116 /***********************************************************************//**
117  * @brief XML constructor
118  *
119  * @param[in] xml XML element.
120  *
121  * Construct spectral nodes model by extracting information from an XML
122  * element. See the read() method for more information about the expected
123  * structure of the XML element.
124  ***************************************************************************/
127 {
128  // Initialise members
129  init_members();
130 
131  // Read information from XML element
132  read(xml);
133 
134  // Return
135  return;
136 }
137 
138 
139 /***********************************************************************//**
140  * @brief Copy constructor
141  *
142  * @param[in] model Spectral nodes model.
143  ***************************************************************************/
145  GModelSpectral(model)
146 {
147  // Initialise members
148  init_members();
149 
150  // Copy members
151  copy_members(model);
152 
153  // Return
154  return;
155 }
156 
157 
158 /***********************************************************************//**
159  * @brief Destructor
160  ***************************************************************************/
162 {
163  // Free members
164  free_members();
165 
166  // Return
167  return;
168 }
169 
170 
171 /*==========================================================================
172  = =
173  = Operators =
174  = =
175  ==========================================================================*/
176 
177 /***********************************************************************//**
178  * @brief Assignment operator
179  *
180  * @param[in] model Spectral nodes model.
181  * @return Spectral nodes model.
182  ***************************************************************************/
184 {
185  // Execute only if object is not identical
186  if (this != &model) {
187 
188  // Copy base class members
189  this->GModelSpectral::operator=(model);
190 
191  // Free members
192  free_members();
193 
194  // Initialise members
195  init_members();
196 
197  // Copy members
198  copy_members(model);
199 
200  } // endif: object was not identical
201 
202  // Return
203  return *this;
204 }
205 
206 
207 /*==========================================================================
208  = =
209  = Public methods =
210  = =
211  ==========================================================================*/
212 
213 /***********************************************************************//**
214  * @brief Clear spectral nodes model
215  ***************************************************************************/
217 {
218  // Free class members (base and derived classes, derived class first)
219  free_members();
221 
222  // Initialise members
224  init_members();
225 
226  // Return
227  return;
228 }
229 
230 
231 /***********************************************************************//**
232  * @brief Clone spectral nodes model
233 ***************************************************************************/
235 {
236  // Clone spectral nodes model
237  return new GModelSpectralNodes(*this);
238 }
239 
240 
241 /***********************************************************************//**
242  * @brief Evaluate function
243  *
244  * @param[in] srcEng True photon energy.
245  * @param[in] srcTime True photon arrival time.
246  * @param[in] gradients Compute gradients?
247  * @return Model value (ph/cm2/s/MeV).
248  *
249  * Computes
250  *
251  * \f[
252  * S_{\rm E}(E | t) =
253  * 10^{(\log {\tt m\_values[i]}) w_{i} +
254  * (\log {\tt m\_values[i+1]}) w_{i+1}}
255  * \f]
256  *
257  * where
258  * - \f${\tt m\_values[i]}\f$ is the intensity of node \f$i\f$,
259  * - \f${\tt m\_values[i+1]}\f$ is the intensity of node \f$i+1\f$,
260  * - \f$w_{i}\f$ is the weighting of node \f$i\f$,
261  * - \f$w_{i+1}\f$ is the weighting of node \f$i+1\f$, and
262  * - \f${\tt m\_energies[i]} <= E <= {\tt m\_energies[i+1]}\f$.
263  *
264  * The weightings \f$w_{i}\f$ and \f$w_{i+1}\f$ are computed by linear
265  * interpolation (in the log-log plane) between the nodes
266  * \f$(\log {\tt m\_energies[i]}, \log{\tt m\_values[i]})\f$
267  * and
268  * \f$(\log {\tt m\_energies[i+1]}, \log{\tt m\_values[i+1]})\f$
269  * to the requested energy \f$\log E\f$.
270  *
271  * If the @p gradients flag is true the method will also compute the partial
272  * derivatives of the model with respect to the parameters using
273  *
274  * \f[
275  * \frac{\delta S_{\rm E}(E | t)}{\delta {\tt m\_values[i]}} =
276  * \frac{S_{\rm E}(E | t) \, w_i}{{\tt m\_values[i]}}
277  * \f]
278  ***************************************************************************/
279 double GModelSpectralNodes::eval(const GEnergy& srcEng,
280  const GTime& srcTime,
281  const bool& gradients) const
282 {
283  // Update evaluation cache
285 
286  // Set interpolation value
288 
289  // Get indices and weights for interpolation
290  int inx_left = m_log_energies.inx_left();
291  int inx_right = m_log_energies.inx_right();
292  double wgt_left = m_log_energies.wgt_left();
293  double wgt_right = m_log_energies.wgt_right();
294 
295  // Interpolate function
296  double exponent = m_log_values[inx_left] * wgt_left +
297  m_log_values[inx_right] * wgt_right;
298 
299  // Compute linear value
300  double value = std::pow(10.0, exponent);
301 
302  // Optionally compute partial derivatives
303  if (gradients) {
304 
305  // Initialise gradients
306  for (int i = 0; i < m_values.size(); ++i) {
307  m_values[i].factor_gradient(0.0);
308  }
309 
310  // Gradient for left node
311  if (m_values[inx_left].is_free() && m_values[inx_left].factor_value() != 0.0) {
312  double grad = value * wgt_left / m_values[inx_left].factor_value();
313  m_values[inx_left].factor_gradient(grad);
314  }
315 
316  // Gradient for right node
317  if (m_values[inx_right].is_free() && m_values[inx_right].factor_value() != 0.0) {
318  double grad = value * wgt_right / m_values[inx_right].factor_value();
319  m_values[inx_right].factor_gradient(grad);
320  }
321 
322  } // endif: computed partial derivatives
323 
324  // Compile option: Check for NaN/Inf
325  #if defined(G_NAN_CHECK)
326  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
327  std::cout << "*** ERROR: GModelSpectralNodes::eval";
328  std::cout << "(srcEng=" << srcEng;
329  std::cout << ", srcTime=" << srcTime << "):";
330  std::cout << " NaN/Inf encountered";
331  std::cout << " (value=" << value;
332  std::cout << ", exponent=" << exponent;
333  std::cout << ")" << std::endl;
334  }
335  #endif
336 
337  // Return
338  return value;
339 }
340 
341 
342 /***********************************************************************//**
343  * @brief Returns model photon flux between [emin, emax] (units: ph/cm2/s)
344  *
345  * @param[in] emin Minimum photon energy.
346  * @param[in] emax Maximum photon energy.
347  * @return Photon flux (ph/cm2/s).
348  *
349  * Computes
350  *
351  * \f[
352  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) dE
353  * \f]
354  *
355  * where
356  * - [@p emin, @p emax] is an energy interval, and
357  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
358  ***************************************************************************/
359 double GModelSpectralNodes::flux(const GEnergy& emin, const GEnergy& emax) const
360 {
361  // Initialise flux
362  double flux = 0.0;
363 
364  // Compute only if integration range is valid
365  if (emin < emax) {
366 
367  // Update flux cache
369 
370  // Get energy range in MeV
371  double e_min = emin.MeV();
372  double e_max = emax.MeV();
373 
374  // Determine left node index for minimum energy
375  m_lin_energies.set_value(e_min);
376  int inx_emin = m_lin_energies.inx_left();
377 
378  // Determine left node index for maximum energy
379  m_lin_energies.set_value(e_max);
380  int inx_emax = m_lin_energies.inx_left();
381 
382  // If both energies are within the same nodes then simply
383  // integrate over the energy interval using the appropriate power
384  // law parameters
385  if (inx_emin == inx_emax) {
386  flux = m_prefactor[inx_emin] *
388  e_max,
389  m_epivot[inx_emin],
390  m_gamma[inx_emin]);
391  }
392 
393  // ... otherwise integrate over the nodes where emin and emax
394  // resides and all the remaining nodes
395  else {
396 
397  // If we are requested to extrapolate beyond first node,
398  // the use the first nodes lower energy and upper energy
399  // boundary for the initial integration.
400  int i_start = (e_min < m_lin_energies[0]) ? inx_emin : inx_emin+1;
401 
402  // Integrate from emin to the node boundary
403  flux = m_prefactor[inx_emin] *
405  m_lin_energies[i_start],
406  m_epivot[inx_emin],
407  m_gamma[inx_emin]);
408 
409  // Integrate over all nodes between
410  for (int i = i_start; i < inx_emax; ++i) {
411  flux += m_flux[i];
412  }
413 
414  // Integrate from node boundary to emax
415  flux += m_prefactor[inx_emax] *
417  e_max,
418  m_epivot[inx_emax],
419  m_gamma[inx_emax]);
420  }
421 
422  } // endif: integration range was valid
423 
424  // Return
425  return flux;
426 }
427 
428 
429 /***********************************************************************//**
430  * @brief Returns model energy flux between [emin, emax] (units: erg/cm2/s)
431  *
432  * @param[in] emin Minimum photon energy.
433  * @param[in] emax Maximum photon energy.
434  * @return Energy flux (erg/cm2/s).
435  *
436  * Computes
437  *
438  * \f[
439  * \int_{\tt emin}^{\tt emax} S_{\rm E}(E | t) E \, dE
440  * \f]
441  *
442  * where
443  * - [@p emin, @p emax] is an energy interval, and
444  * - \f$S_{\rm E}(E | t)\f$ is the spectral model (ph/cm2/s/MeV).
445  ***************************************************************************/
446 double GModelSpectralNodes::eflux(const GEnergy& emin, const GEnergy& emax) const
447 {
448  // Initialise flux
449  double eflux = 0.0;
450 
451  // Compute only if integration range is valid
452  if (emin < emax) {
453 
454  // Update flux cache
456 
457  // Get energy range in MeV
458  double e_min = emin.MeV();
459  double e_max = emax.MeV();
460 
461  // Determine left node index for minimum energy
462  m_lin_energies.set_value(e_min);
463  int inx_emin = m_lin_energies.inx_left();
464 
465  // Determine left node index for maximum energy
466  m_lin_energies.set_value(e_max);
467  int inx_emax = m_lin_energies.inx_left();
468 
469  // If both energies are within the same nodes then simply
470  // integrate over the energy interval using the appropriate power
471  // law parameters
472  if (inx_emin == inx_emax) {
473  eflux = m_prefactor[inx_emin] *
475  e_max,
476  m_epivot[inx_emin],
477  m_gamma[inx_emin]) *
479  }
480 
481  // ... otherwise integrate over the nodes where emin and emax
482  // resides and all the remaining nodes
483  else {
484 
485  // If we are requested to extrapolate beyond first node,
486  // the use the first nodes lower energy and upper energy
487  // boundary for the initial integration.
488  int i_start = (e_min < m_lin_energies[0]) ? inx_emin : inx_emin+1;
489 
490  // Integrate from emin to the node boundary
491  eflux = m_prefactor[inx_emin] *
493  m_lin_energies[i_start],
494  m_epivot[inx_emin],
495  m_gamma[inx_emin]) *
497 
498  // Integrate over all nodes between
499  for (int i = i_start; i < inx_emax; ++i) {
500  eflux += m_eflux[i];
501  }
502 
503  // Integrate from node boundary to emax
504  eflux += m_prefactor[inx_emax] *
506  e_max,
507  m_epivot[inx_emax],
508  m_gamma[inx_emax]) *
510 
511  } // endelse: emin and emax not between same nodes
512 
513  } // endif: integration range was valid
514 
515  // Return flux
516  return eflux;
517 }
518 
519 
520 /***********************************************************************//**
521  * @brief Returns MC energy between [emin, emax]
522  *
523  * @param[in] emin Minimum photon energy.
524  * @param[in] emax Maximum photon energy.
525  * @param[in] time True photon arrival time.
526  * @param[in,out] ran Random number generator.
527  * @return Energy.
528  *
529  * @exception GException::invalid_return_value
530  * No valid Monte Carlo cache
531  *
532  * Returns Monte Carlo energy by randomly drawing from node function.
533  ***************************************************************************/
535  const GEnergy& emax,
536  const GTime& time,
537  GRan& ran) const
538 {
539  // Check energy interval
541 
542  // Allocate energy
543  GEnergy energy;
544 
545  // Update cache
546  mc_update(emin, emax);
547 
548  // Determine in which bin we reside. If this fails then thrown an
549  // exception
550  int inx = 0;
551  if (m_mc_cum.size() > 1) {
552  double u = ran.uniform();
553  for (inx = m_mc_cum.size()-1; inx > 0; --inx) {
554  if (m_mc_cum[inx-1] <= u) {
555  break;
556  }
557  }
558  }
559  else if (m_mc_cum.size() == 0) {
560  std::string msg = "No valid nodes found for energy interval ["+
561  emin.print()+","+emax.print()+"]. Either restrict "
562  "the energy range to the one covered by the "
563  "spectral nodes or extend the spectral nodes "
564  "in energy.";
566  }
567 
568  // Get random energy for specific bin
569  if (m_mc_exp[inx] != 0.0) {
570  double e_min = m_mc_min[inx];
571  double e_max = m_mc_max[inx];
572  double u = ran.uniform();
573  double eng = (u > 0.0)
574  ? std::exp(std::log(u * (e_max - e_min) + e_min) / m_mc_exp[inx])
575  : 0.0;
576  energy.MeV(eng);
577  }
578  else {
579  double e_min = m_mc_min[inx];
580  double e_max = m_mc_max[inx];
581  double u = ran.uniform();
582  double eng = std::exp(u * (e_max - e_min) + e_min);
583  energy.MeV(eng);
584  }
585 
586  // Return energy
587  return energy;
588 }
589 
590 
591 /***********************************************************************//**
592  * @brief Read model from XML element
593  *
594  * @param[in] xml XML element.
595  *
596  * @exception GException::invalid_value
597  * Not enough nodes in node function.
598  * @exception GException::invalid_value
599  * Energy or intensity are not positive.
600  *
601  * Reads the spectral information from an XML element. The format of the XML
602  * elements is
603  *
604  * <spectrum type="NodeFunction">
605  * <node>
606  * <parameter name="Energy" ../>
607  * <parameter name="Intensity" ../>
608  * </node>
609  * ...
610  * <node>
611  * <parameter name="Energy" ../>
612  * <parameter name="Intensity" ../>
613  * </node>
614  * </spectrum>
615  *
616  * @todo Check that nodes are ordered
617  * @todo Check that energy boundaries are not overlapping
618  * @todo Check whether we require at least one or two nodes
619  ***************************************************************************/
621 {
622  // Free space for nodes
623  m_energies.clear();
624  m_values.clear();
625 
626  // Get number of nodes from XML file
627  int nodes = xml.elements("node");
628 
629  // Throw an error if there not at least two nodes
630  if (nodes < 1) {
631  std::string msg = "Node function model contains "+gammalib::str(nodes)+
632  " but requires at least one node. Please specify at "
633  "least one node.";
634  throw GException::invalid_value(G_READ, msg);
635  }
636 
637  // Loop over all nodes
638  for (int i = 0; i < nodes; ++i) {
639 
640  // Allocate node parameters
643 
644  // Get node
645  const GXmlElement* node = xml.element("node", i);
646 
647  // Get parameters
648  const GXmlElement* epar = gammalib::xml_get_par(G_READ, *node, "Energy");
649  const GXmlElement* ipar = gammalib::xml_get_par(G_READ, *node, "Intensity");
650 
651  // Read parameters
652  energy.read(*epar);
653  intensity.read(*ipar);
654 
655  // Throw an exception if either energy or intensity is not positive
656  if (energy.value() <= 0.0) {
657  std::string msg = "Non positive energy "+
658  gammalib::str(energy.value())+" MeV encountered "
659  "in model definition XML file. Please specify "
660  "only nodes with positive energy.";
661  throw GException::invalid_value(G_READ, msg);
662  }
663  if (intensity.value() <= 0.0) {
664  std::string msg = "Non positive intensity "+
665  gammalib::str(intensity.value())+" ph/cm2/s/MeV "
666  "encountered in model definition XML file. "
667  "Please specify only nodes with positive "
668  "intensity.";
669  throw GException::invalid_value(G_READ, msg);
670  }
671 
672  // Set parameter names
673  std::string energy_name = "Energy"+gammalib::str(i);
674  std::string intensity_name = "Intensity"+gammalib::str(i);
675 
676  // Set energy attributes
677  energy.name(energy_name);
678  energy.unit("MeV");
679  energy.has_grad(false);
680 
681  // Set intensity attributes
682  intensity.name(intensity_name);
683  intensity.unit("ph/cm2/s/MeV");
684  intensity.has_grad(true);
685 
686  // Push node parameters on list
687  m_energies.push_back(energy);
688  m_values.push_back(intensity);
689 
690  } // endfor: looped over nodes
691 
692  // Update parameter mapping
693  update_pars();
694 
695  // Set pre-computation cache
696  set_cache();
697 
698  // Return
699  return;
700 }
701 
702 
703 /***********************************************************************//**
704  * @brief Write model into XML element
705  *
706  * @param[in] xml XML element into which model information is written.
707  *
708  * @exception GException::invalid_value
709  * Invalid number of nodes found in XML element.
710  *
711  * Writes the spectral information into an XML element. The format of the XML
712  * element is
713  *
714  * <spectrum type="NodeFunction">
715  * <node>
716  * <parameter name="Energy" ../>
717  * <parameter name="Intensity" ../>
718  * </node>
719  * ...
720  * <node>
721  * <parameter name="Energy" ../>
722  * <parameter name="Intensity" ../>
723  * </node>
724  * </spectrum>
725  ***************************************************************************/
727 {
728  // Determine number of nodes
729  int nodes = m_energies.size();
730 
731  // Verify model type
733 
734  // If XML element has 0 nodes then append nodes
735  if (xml.elements() == 0) {
736  for (int i = 0; i < nodes; ++i) {
737  xml.append(GXmlElement("node"));
738  }
739  }
740 
741  // Verify that XML element has the required number of nodes
742  if (xml.elements() != nodes) {
743  std::string msg = "Number of "+gammalib::str(xml.elements())+
744  " child elements in XML file does not correspond "
745  "to expected number of "+gammalib::str(nodes)+
746  " elements. Please verify the XML format.";
748  }
749  int npars = xml.elements("node");
750  if (npars != nodes) {
751  std::string msg = "Number of "+gammalib::str(npars)+" \"node\" "
752  "child elements in XML file does not correspond to "
753  "expected number of "+gammalib::str(nodes)+
754  " elements. Please verify the XML format.";
756  }
757 
758  // Loop over all nodes
759  for (int i = 0; i < nodes; ++i) {
760 
761  // Get node
762  GXmlElement* node = xml.element("node", i);
763 
764  // Get copy of parameters so that we can change their names
767 
768  // Set names since we appended for internal usage the indices to the
769  // parameter names, and we want to get rid of them for writing the
770  // model into the XML element
771  energy.name("Energy");
772  intensity.name("Intensity");
773 
774  // Get XML parameters
775  GXmlElement* eng = gammalib::xml_need_par(G_WRITE, *node, energy.name());
776  GXmlElement* val = gammalib::xml_need_par(G_WRITE, *node, intensity.name());
777 
778  // Write parameters
779  energy.write(*eng);
780  intensity.write(*val);
781 
782  } // endfor: looped over nodes
783 
784  // Return
785  return;
786 }
787 
788 
789 /***********************************************************************//**
790  * @brief Append node
791  *
792  * @param[in] energy Node energy.
793  * @param[in] intensity Node intensity.
794  *
795  * @exception GException::invalid_argument
796  * Non-positive energy or intensity specified.
797  *
798  * Appends one node to the node function. By default the energy of the node
799  * is fixed while the intensity of the node is free.
800  ***************************************************************************/
802  const double& intensity)
803 {
804  // Throw an exception if either energy or intensity is not positive
805  if (energy.MeV() <= 0.0) {
806  std::string msg = "Non-positive energy "+energy.print()+" not allowed. "
807  "Please specify only positive energies.";
809  }
810  if (intensity <= 0.0) {
811  std::string msg = "Non-positive intensity "+
812  gammalib::str(intensity)+" ph/cm2/s/MeV "
813  "not allowed. Please specify only positive "
814  "intensities.";
816  }
817 
818  // Allocate node parameters
819  GModelPar e_par;
820  GModelPar i_par;
821 
822  // Set energy attributes
823  e_par.name("Energy");
824  e_par.value(energy.MeV());
825  e_par.unit("MeV");
826  e_par.has_grad(false);
827  e_par.fix();
828 
829  // Set intensity attributes
830  i_par.name("Intensity");
831  i_par.value(intensity);
832  i_par.unit("ph/cm2/s/MeV");
833  i_par.has_grad(true);
834  i_par.free();
835 
836  // Append to nodes
837  m_energies.push_back(e_par);
838  m_values.push_back(i_par);
839 
840  // Update parameter mapping
841  update_pars();
842 
843  // Set pre-computation cache
844  set_cache();
845 
846  // Return
847  return;
848 }
849 
850 
851 /***********************************************************************//**
852  * @brief Insert node
853  *
854  * @param[in] index Node index [0,...,nodes()[.
855  * @param[in] energy Node energy.
856  * @param[in] intensity Node intensity.
857  *
858  * @exception GException::out_of_range
859  * Node index is out of range.
860  * @exception GException::invalid_argument
861  * Non-positive energy or intensity specified.
862  *
863  * Inserts a node into the node function before the node with the specified
864  * @p index. By default the energy of the node is fixed while the intensity
865  * of the node is free.
866  ***************************************************************************/
867 void GModelSpectralNodes::insert(const int& index,
868  const GEnergy& energy,
869  const double& intensity)
870 {
871  // Raise exception if index is outside boundary
872  #if defined(G_RANGE_CHECK)
873  if (m_energies.empty()) {
874  if (index > 0) {
875  throw GException::out_of_range(G_INSERT, "Node index",
876  index, nodes());
877  }
878  }
879  else {
880  if (index < 0 || index >= nodes()) {
881  throw GException::out_of_range(G_INSERT, "Node index",
882  index, nodes());
883  }
884  }
885  #endif
886 
887  // Throw an exception if either energy or intensity is not positive
888  if (energy.MeV() <= 0.0) {
889  std::string msg = "Non-positive energy "+energy.print()+" not allowed. "
890  "Please specify only positive energies.";
892  }
893  if (intensity <= 0.0) {
894  std::string msg = "Non-positive intensity "+
895  gammalib::str(intensity)+" ph/cm2/s/MeV "
896  "not allowed. Please specify only positive "
897  "intensities.";
899  }
900 
901  // Allocate node parameters
902  GModelPar e_par;
903  GModelPar i_par;
904 
905  // Set energy attributes
906  e_par.name("Energy");
907  e_par.value(energy.MeV());
908  e_par.unit("MeV");
909  e_par.has_grad(false);
910  e_par.fix();
911 
912  // Set intensity attributes
913  i_par.name("Intensity");
914  i_par.value(intensity);
915  i_par.unit("ph/cm2/s/MeV");
916  i_par.has_grad(true);
917  i_par.free();
918 
919  // Insert node
920  m_energies.insert(m_energies.begin()+index, e_par);
921  m_values.insert(m_values.begin()+index, i_par);
922 
923  // Update parameter mapping
924  update_pars();
925 
926  // Set pre-computation cache
927  set_cache();
928 
929  // Return
930  return;
931 }
932 
933 
934 /***********************************************************************//**
935  * @brief Remove node
936  *
937  * @param[in] index Node index [0,...,nodes()[.
938  *
939  * @exception GException::out_of_range
940  * Node index is out of range.
941  *
942  * Removes node of specified @p index from the node function.
943  ***************************************************************************/
944 void GModelSpectralNodes::remove(const int& index)
945 {
946  // Raise exception if index is outside boundary
947  #if defined(G_RANGE_CHECK)
948  if (index < 0 || index >= nodes()) {
949  throw GException::out_of_range(G_REMOVE, "Node index",
950  index, nodes());
951  }
952  #endif
953 
954  // Erase energy and intensity
955  m_energies.erase(m_energies.begin() + index);
956  m_values.erase(m_values.begin() + index);
957 
958  // Update parameter mapping
959  update_pars();
960 
961  // Set pre-computation cache
962  set_cache();
963 
964  // Return
965  return;
966 }
967 
968 
969 /***********************************************************************//**
970  * @brief Reserve space for nodes
971  *
972  * @param[in] num Number of reseved nodes.
973  *
974  * Reserves space for @p num nodes.
975  ***************************************************************************/
976 void GModelSpectralNodes::reserve(const int& num)
977 {
978  // Reserve space
979  m_energies.reserve(num);
980  m_values.reserve(num);
981 
982  // Return
983  return;
984 }
985 
986 
987 /***********************************************************************//**
988  * @brief Append nodes from node function
989  *
990  * @param[in] nodes Node function.
991  *
992  * Appends all nodes from a node function to current object.
993  ***************************************************************************/
995 {
996  // Get number of nodes in node function. Note that we extract the size
997  // first to avoid an endless loop that arises when a container is
998  // appended to itself.
999  int num = nodes.nodes();
1000 
1001  // Continue only if node function is not empty
1002  if (num > 0) {
1003 
1004  // Reserve enough space
1005  reserve(this->nodes() + num);
1006 
1007  // Loop over all nodes and append them to the node function
1008  for (int i = 0; i < num; ++i) {
1009  m_energies.push_back(nodes.m_energies[i]);
1010  m_values.push_back(nodes.m_values[i]);
1011  }
1012 
1013  // Update parameter mapping
1014  update_pars();
1015 
1016  // Set pre-computation cache
1017  set_cache();
1018 
1019  } // endif: node function was not empty
1020 
1021  // Return
1022  return;
1023 }
1024 
1025 
1026 /***********************************************************************//**
1027  * @brief Return node energy
1028  *
1029  * @param[in] index Node index [0,...,nodes()[.
1030  * @return Energy of node @p index.
1031  *
1032  * @exception GException::out_of_range
1033  * Index is out of range.
1034  *
1035  * Returns the energy of node @p index.
1036  ***************************************************************************/
1037 GEnergy GModelSpectralNodes::energy(const int& index) const
1038 {
1039  // Raise an exception if index is out of range
1040  #if defined(G_RANGE_CHECK)
1041  if (index < 0 || index >= nodes()) {
1042  throw GException::out_of_range(G_ENERGY_GET, "Node index",
1043  index, nodes());
1044  }
1045  #endif
1046 
1047  // Retrieve energy
1048  GEnergy energy;
1049  energy.MeV(m_energies[index].value());
1050 
1051  // Return energy
1052  return energy;
1053 }
1054 
1055 
1056 /***********************************************************************//**
1057  * @brief Set node energy
1058  *
1059  * @param[in] index Node index [0,...,nodes()[.
1060  * @param[in] energy Node energy.
1061  *
1062  * @exception GException::out_of_range
1063  * Index is out of range.
1064  * @exception GException::invalid_argument
1065  * Non-positive energy specified.
1066  *
1067  * Sets the energy of node @p index.
1068  ***************************************************************************/
1069 void GModelSpectralNodes::energy(const int& index, const GEnergy& energy)
1070 {
1071  // Raise an exception if index is out of range
1072  #if defined(G_RANGE_CHECK)
1073  if (index < 0 || index >= nodes()) {
1074  throw GException::out_of_range(G_ENERGY_SET, "Node index",
1075  index, nodes());
1076  }
1077  #endif
1078 
1079  // Throw an exception if energy is not positive
1080  if (energy.MeV() <= 0.0) {
1081  std::string msg = "Non-positive energy "+energy.print()+" not allowed. "
1082  "Please specify only positive energies.";
1084  }
1085 
1086  // Set energy
1087  m_energies[index].value(energy.MeV());
1088 
1089  // Set pre-computation cache
1090  set_cache();
1091 
1092  // Return
1093  return;
1094 }
1095 
1096 
1097 /***********************************************************************//**
1098  * @brief Return node intensity
1099  *
1100  * @param[in] index Node index [0,...,nodes()[.
1101  * @return Intensity of node @p index.
1102  *
1103  * @exception GException::out_of_range
1104  * Index is out of range.
1105  *
1106  * Returns the intensity of node @p index.
1107  ***************************************************************************/
1108 double GModelSpectralNodes::intensity(const int& index) const
1109 {
1110  // Raise an exception if index is out of range
1111  #if defined(G_RANGE_CHECK)
1112  if (index < 0 || index >= nodes()) {
1113  throw GException::out_of_range(G_INTENSITY_GET, "Node index",
1114  index, nodes());
1115  }
1116  #endif
1117 
1118  // Return intensity
1119  return (m_values[index].value());
1120 }
1121 
1122 
1123 /***********************************************************************//**
1124  * @brief Return intensity error of node
1125  *
1126  * @param[in] index Node index [0,...,nodes()[.
1127  * @return Intensity error of node @p index.
1128  *
1129  * @exception GException::out_of_range
1130  * Index is out of range.
1131  *
1132  * Returns the intensity error of node @p index.
1133  ***************************************************************************/
1134 double GModelSpectralNodes::error(const int& index) const
1135 {
1136  // Raise an exception if index is out of range
1137  #if defined(G_RANGE_CHECK)
1138  if (index < 0 || index >= nodes()) {
1139  throw GException::out_of_range(G_ERROR_GET, "Node index",
1140  index, nodes());
1141  }
1142  #endif
1143 
1144  // Return intensity error
1145  return (m_values[index].error());
1146 }
1147 
1148 
1149 /***********************************************************************//**
1150  * @brief Set node intensity
1151  *
1152  * @param[in] index Node index [0,...,nodes()[.
1153  * @param[in] intensity Node Intensity.
1154  *
1155  * @exception GException::out_of_range
1156  * Index is out of range.
1157  * @exception GException::invalid_argument
1158  * Non-positive intensity specified.
1159  *
1160  * Set the intensity of node @p index.
1161  ***************************************************************************/
1162 void GModelSpectralNodes::intensity(const int& index, const double& intensity)
1163 {
1164  // Raise an exception if index is out of range
1165  #if defined(G_RANGE_CHECK)
1166  if (index < 0 || index >= nodes()) {
1167  throw GException::out_of_range(G_INTENSITY_SET, "Node index",
1168  index, nodes());
1169  }
1170  #endif
1171 
1172  // Throw an exception if either energy or intensity is not positive
1173  if (intensity <= 0.0) {
1174  std::string msg = "Non-positive intensity "+
1175  gammalib::str(intensity)+" ph/cm2/s/MeV "
1176  "not allowed. Please specify only positive "
1177  "intensities.";
1179  }
1180 
1181  // Set intensity
1182  m_values[index].value(intensity);
1183 
1184  // Set pre-computation cache
1185  set_cache();
1186 
1187  // Return
1188  return;
1189 }
1190 
1191 
1192 /***********************************************************************//**
1193  * @brief Print node function information
1194  *
1195  * @param[in] chatter Chattiness.
1196  * @return String containing node function information.
1197  ***************************************************************************/
1198 std::string GModelSpectralNodes::print(const GChatter& chatter) const
1199 {
1200  // Initialise result string
1201  std::string result;
1202 
1203  // Continue only if chatter is not silent
1204  if (chatter != SILENT) {
1205 
1206  // Append header
1207  result.append("=== GModelSpectralNodes ===");
1208 
1209  // Append information
1210  result.append("\n"+gammalib::parformat("Number of nodes"));
1211  result.append(gammalib::str(m_energies.size()));
1212  result.append("\n"+gammalib::parformat("Number of parameters"));
1213  result.append(gammalib::str(size()));
1214  for (int i = 0; i < size(); ++i) {
1215  result.append("\n"+m_pars[i]->print(chatter));
1216  }
1217 
1218  // VERBOSE: Append information about power laws between node
1219  if (chatter == VERBOSE) {
1220  for (int i = 0; i < m_prefactor.size(); ++i) {
1221  result.append("\n"+gammalib::parformat("Interval "+gammalib::str(i+1)));
1222  result.append("Epivot="+gammalib::str(m_epivot[i]));
1223  result.append(" Prefactor="+gammalib::str(m_prefactor[i]));
1224  result.append(" Gamma="+gammalib::str(m_gamma[i]));
1225  result.append(" Flux="+gammalib::str(m_flux[i]));
1226  result.append(" EFlux="+gammalib::str(m_eflux[i]));
1227  }
1228  }
1229 
1230  } // endif: chatter was not silent
1231 
1232  // Return result
1233  return result;
1234 }
1235 
1236 
1237 /*==========================================================================
1238  = =
1239  = Private methods =
1240  = =
1241  ==========================================================================*/
1242 
1243 /***********************************************************************//**
1244  * @brief Initialise class members
1245  ***************************************************************************/
1247 {
1248  // Initialise node energies and values
1249  m_energies.clear();
1250  m_values.clear();
1251 
1252  // Initialise evaluation cache
1253  m_old_energies.clear();
1254  m_old_values.clear();
1256  m_log_values.clear();
1257 
1258  // Initialise flux computation cache
1260  m_lin_values.clear();
1261  m_prefactor.clear();
1262  m_gamma.clear();
1263  m_epivot.clear();
1264  m_flux.clear();
1265  m_eflux.clear();
1266 
1267  // Initialise MC cache
1268  m_mc_emin.clear();
1269  m_mc_emax.clear();
1270  m_mc_cum.clear();
1271  m_mc_min.clear();
1272  m_mc_max.clear();
1273  m_mc_exp.clear();
1274 
1275  // Update parameter mapping
1276  update_pars();
1277 
1278  // Return
1279  return;
1280 }
1281 
1282 
1283 /***********************************************************************//**
1284  * @brief Copy class members
1285  *
1286  * @param[in] model Spectral function model.
1287  ***************************************************************************/
1289 {
1290  // Copy node energies and values
1291  m_energies = model.m_energies;
1292  m_values = model.m_values;
1293 
1294  // Copy evaluation cache
1296  m_old_values = model.m_old_values;
1298  m_log_values = model.m_log_values;
1299 
1300  // Copy flux computation cache
1302  m_lin_values = model.m_lin_values;
1303  m_prefactor = model.m_prefactor;
1304  m_gamma = model.m_gamma;
1305  m_epivot = model.m_epivot;
1306  m_flux = model.m_flux;
1307  m_eflux = model.m_eflux;
1308 
1309  // Copy MC cache
1310  m_mc_emin = model.m_mc_emin;
1311  m_mc_emax = model.m_mc_emax;
1312  m_mc_cum = model.m_mc_cum;
1313  m_mc_min = model.m_mc_min;
1314  m_mc_max = model.m_mc_max;
1315  m_mc_exp = model.m_mc_exp;
1316 
1317  // Update parameter mapping
1318  update_pars();
1319 
1320  // Return
1321  return;
1322 }
1323 
1324 
1325 /***********************************************************************//**
1326  * @brief Delete class members
1327  ***************************************************************************/
1329 {
1330  // Return
1331  return;
1332 }
1333 
1334 
1335 /***********************************************************************//**
1336  * @brief Update parameter mapping
1337  *
1338  * Sets the parameter pointers in the m_pars array, enabling iterating over
1339  * all model parameters. This method needs to be called after changing the
1340  * number of nodes in the spectral model. The method needs not to be called
1341  * after value update.
1342  ***************************************************************************/
1344 {
1345  // Clear parameter pointer(s)
1346  m_pars.clear();
1347 
1348  // Get number of nodes
1349  int nodes = m_energies.size();
1350 
1351  // Set parameter pointers for all nodes
1352  for (int i = 0; i < nodes; ++i) {
1353 
1354  // Set parameter names
1355  std::string energy_name = "Energy"+gammalib::str(i);
1356  std::string intensity_name = "Intensity"+gammalib::str(i);
1357 
1358  // Set energy attributes
1359  m_energies[i].name(energy_name);
1360  m_energies[i].unit("MeV");
1361  m_energies[i].has_grad(false);
1362 
1363  // Set intensity attributes
1364  m_values[i].name(intensity_name);
1365  m_values[i].unit("ph/cm2/s/MeV");
1366  m_values[i].has_grad(true);
1367 
1368  // Set pointer
1369  m_pars.push_back(&(m_energies[i]));
1370  m_pars.push_back(&(m_values[i]));
1371 
1372  } // endfor: looped over nodes
1373 
1374  // Return
1375  return;
1376 }
1377 
1378 
1379 /***********************************************************************//**
1380  * @brief Set pre-computation cache
1381  ***************************************************************************/
1383 {
1384  // Set evaluation cache
1385  set_eval_cache();
1386 
1387  // Set flux computation cache
1388  set_flux_cache();
1389 
1390  // Return
1391  return;
1392 }
1393 
1394 
1395 /***********************************************************************//**
1396  * @brief Set evaluation cache
1397  *
1398  * The evaluation cache contains pre-computed results that are needed for
1399  * fast function evaluation.
1400  *
1401  * @todo Check that all energies and intensities are > 0
1402  ***************************************************************************/
1404 {
1405  // Clear any existing values
1406  m_old_energies.clear();
1407  m_old_values.clear();
1409  m_log_values.clear();
1410 
1411  // Compute log10 of energies and intensities
1412  for (int i = 0; i < m_energies.size(); ++i) {
1413 
1414  // Set log10(energy)
1415  double log10energy = std::log10(m_energies[i].value());
1416  m_old_energies.push_back(log10energy);
1417  m_log_energies.append(log10energy);
1418 
1419  // Set log10(intensity)
1420  double log10value = std::log10(m_values[i].value());
1421  m_old_values.push_back(log10value);
1422  m_log_values.push_back(log10value);
1423 
1424  }
1425 
1426  // Return
1427  return;
1428 }
1429 
1430 
1431 /***********************************************************************//**
1432  * @brief Set flux computation cache
1433  *
1434  * Pre-computes some values that are needed for flux computation.
1435  *
1436  * @todo Handle special case emin=emax and fmin=fmax
1437  ***************************************************************************/
1439 {
1440  // Clear any existing values
1442  m_lin_values.clear();
1443  m_prefactor.clear();
1444  m_gamma.clear();
1445  m_epivot.clear();
1446  m_flux.clear();
1447  m_eflux.clear();
1448 
1449  // Store linear energies and values
1450  for (int i = 0; i < m_energies.size(); ++i) {
1451  m_lin_energies.append(m_energies[i].value());
1452  m_lin_values.push_back(m_values[i].value());
1453  }
1454 
1455  // Loop over all nodes-1
1456  for (int i = 0; i < m_energies.size()-1; ++i) {
1457 
1458  // Get energies and function values
1459  double emin = m_lin_energies[i];
1460  double emax = m_lin_energies[i+1];
1461  double fmin = m_values[i].value();
1462  double fmax = m_values[i+1].value();
1463 
1464  // Compute pivot energy (MeV). We use here the geometric mean of
1465  // the node boundaries.
1466  double epivot = std::sqrt(emin*emax);
1467 
1468  // Compute spectral index
1469  double gamma = std::log(fmin/fmax) / std::log(emin/emax);
1470 
1471  // Compute power law normalisation
1472  double prefactor = fmin / std::pow(emin/epivot, gamma);
1473 
1474  // Compute photon flux between nodes
1475  double flux = prefactor *
1476  gammalib::plaw_photon_flux(emin, emax, epivot, gamma);
1477 
1478  // Compute energy flux between nodes
1479  double eflux = prefactor *
1480  gammalib::plaw_energy_flux(emin, emax, epivot, gamma);
1481 
1482  // Convert energy flux from MeV/cm2/s to erg/cm2/s
1483  eflux *= gammalib::MeV2erg;
1484 
1485  // Push back values on pre-computation cache
1486  m_prefactor.push_back(prefactor);
1487  m_gamma.push_back(gamma);
1488  m_epivot.push_back(epivot);
1489  m_flux.push_back(flux);
1490  m_eflux.push_back(eflux);
1491 
1492  } // endfor: looped over all nodes
1493 
1494  // Return
1495  return;
1496 }
1497 
1498 
1499 /***********************************************************************//**
1500  * @brief Update evaluation cache
1501  *
1502  * Updates the evaluation cache by computing only results for values that
1503  * changed.
1504  *
1505  * @todo Check that all energies and intensities are > 0
1506  ***************************************************************************/
1508 {
1509  // Update energies
1510  for (int i = 0; i < m_energies.size(); ++i) {
1511  double energy = m_energies[i].value();
1512  if (energy != m_old_energies[i]) {
1513  m_log_energies[i] = std::log10(energy);
1514  m_old_energies[i] = energy;
1515  }
1516  }
1517 
1518  // Update intensities
1519  for (int i = 0; i < m_values.size(); ++i) {
1520  double value = m_values[i].value();
1521  if (value != m_old_values[i]) {
1522  m_log_values[i] = std::log10(value);
1523  m_old_values[i] = value;
1524  }
1525  }
1526 
1527  // Return
1528  return;
1529 }
1530 
1531 
1532 /***********************************************************************//**
1533  * @brief Update flux computation cache
1534  *
1535  * Updates the flux computation cache if either the energy boundaries or the
1536  * intensity values have changed.
1537  ***************************************************************************/
1539 {
1540  // Loop over all nodes. If energy or value of a node has changed then
1541  // set cache
1542  for (int i = 0; i < m_energies.size(); ++i) {
1543  if ((m_lin_energies[i] != m_energies[i].value()) ||
1544  (m_lin_values[i] != m_values[i].value())) {
1545  set_flux_cache();
1546  break;
1547  }
1548  }
1549 
1550  // Return
1551  return;
1552 }
1553 
1554 
1555 /***********************************************************************//**
1556  * @brief Set MC pre-computation cache
1557  *
1558  * @param[in] emin Minimum energy.
1559  * @param[in] emax Maximum energy.
1560  *
1561  * This method sets up an array of indices and the cumulative distribution
1562  * function needed for MC simulations.
1563  ***************************************************************************/
1564 void GModelSpectralNodes::mc_update(const GEnergy& emin, const GEnergy& emax) const
1565 {
1566  // Check if we need to update the cache
1567  if (emin != m_mc_emin || emax != m_mc_emax) {
1568 
1569  // Store new energy interval
1570  m_mc_emin = emin;
1571  m_mc_emax = emax;
1572 
1573  // Initialise cache
1574  m_mc_cum.clear();
1575  m_mc_min.clear();
1576  m_mc_max.clear();
1577  m_mc_exp.clear();
1578 
1579  // Get energy range in MeV
1580  double e_min = emin.MeV();
1581  double e_max = emax.MeV();
1582 
1583  // Continue only if e_max > e_min
1584  if (e_max > e_min) {
1585 
1586  // Allocate flux
1587  double flux;
1588 
1589  // Determine left node index for minimum energy
1590  m_lin_energies.set_value(e_min);
1591  int inx_emin = m_lin_energies.inx_left();
1592 
1593  // Determine left node index for maximum energy
1594  m_lin_energies.set_value(e_max);
1595  int inx_emax = m_lin_energies.inx_left();
1596 
1597  // If both energies are within the same node then just
1598  // add this one node on the stack
1599  if (inx_emin == inx_emax) {
1600  flux = m_prefactor[inx_emin] *
1602  e_max,
1603  m_epivot[inx_emin],
1604  m_gamma[inx_emin]);
1605  m_mc_cum.push_back(flux);
1606  m_mc_min.push_back(e_min);
1607  m_mc_max.push_back(e_max);
1608  m_mc_exp.push_back(m_gamma[inx_emin]);
1609  }
1610 
1611  // ... otherwise integrate over the nodes where emin and emax
1612  // resides and all the remaining nodes
1613  else {
1614 
1615  // If we are requested to extrapolate beyond first node,
1616  // the use the first nodes lower energy and upper energy
1617  // boundary for the initial integration.
1618  int i_start = (e_min < m_lin_energies[0]) ? inx_emin : inx_emin+1;
1619 
1620  // Add emin to the node boundary
1621  flux = m_prefactor[inx_emin] *
1623  m_lin_energies[i_start],
1624  m_epivot[inx_emin],
1625  m_gamma[inx_emin]);
1626  m_mc_cum.push_back(flux);
1627  m_mc_min.push_back(e_min);
1628  m_mc_max.push_back(m_lin_energies[i_start]);
1629  m_mc_exp.push_back(m_gamma[inx_emin]);
1630 
1631  // Add all nodes between
1632  for (int i = i_start; i < inx_emax; ++i) {
1633  flux = m_flux[i];
1634  m_mc_cum.push_back(flux);
1635  m_mc_min.push_back(m_lin_energies[i]);
1636  m_mc_max.push_back(m_lin_energies[i+1]);
1637  m_mc_exp.push_back(m_gamma[i]);
1638  }
1639 
1640  // Add node boundary to emax
1641  flux = m_prefactor[inx_emax] *
1643  e_max,
1644  m_epivot[inx_emax],
1645  m_gamma[inx_emax]);
1646  m_mc_cum.push_back(flux);
1647  m_mc_min.push_back(m_lin_energies[inx_emax]);
1648  m_mc_max.push_back(e_max);
1649  m_mc_exp.push_back(m_gamma[inx_emax]);
1650 
1651  } // endelse: emin and emax not between same nodes
1652 
1653  // Build cumulative distribution
1654  for (int i = 1; i < m_mc_cum.size(); ++i) {
1655  m_mc_cum[i] += m_mc_cum[i-1];
1656  }
1657  double norm = m_mc_cum[m_mc_cum.size()-1];
1658  if (norm > 0.0) {
1659  for (int i = 0; i < m_mc_cum.size(); ++i) {
1660  m_mc_cum[i] /= norm;
1661  }
1662  }
1663 
1664  // Set MC values
1665  for (int i = 0; i < m_mc_cum.size(); ++i) {
1666 
1667  // Compute exponent
1668  double exponent = m_mc_exp[i] + 1.0;
1669 
1670  // Exponent dependend computation
1671  if (std::abs(exponent) > 1.0e-11) {
1672 
1673  // If the exponent is too small then use lower energy
1674  // boundary
1675  if (exponent < -50.0) {
1676  m_mc_exp[i] = 0.0;
1677  m_mc_min[i] = std::log(m_mc_min[i]);
1678  m_mc_max[i] = m_mc_min[i];
1679  }
1680 
1681  // ... otherwise if exponent is too large then use
1682  // upper energy boundary
1683  else if (exponent > +50.0) {
1684  m_mc_exp[i] = 0.0;
1685  m_mc_min[i] = std::log(m_mc_max[i]);
1686  m_mc_max[i] = m_mc_min[i];
1687  }
1688 
1689  // ... otherwise use transformation formula
1690  else {
1691  m_mc_exp[i] = exponent;
1692  m_mc_min[i] = std::pow(m_mc_min[i], exponent);
1693  m_mc_max[i] = std::pow(m_mc_max[i], exponent);
1694  }
1695  }
1696  else {
1697  m_mc_exp[i] = 0.0;
1698  m_mc_min[i] = std::log(m_mc_min[i]);
1699  m_mc_max[i] = std::log(m_mc_max[i]);
1700  }
1701 
1702  } // endfor: set MC values
1703 
1704  } // endif: e_max > e_min
1705 
1706  } // endif: Update was required
1707 
1708  // Return
1709  return;
1710 }
std::vector< double > m_flux
Photon fluxes.
void update_flux_cache(void) const
Update flux computation cache.
virtual GEnergy mc(const GEnergy &emin, const GEnergy &emax, const GTime &time, GRan &ran) const
Returns MC energy between [emin, emax].
std::vector< double > m_eflux
Energy fluxes.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
virtual GModelSpectralNodes & operator=(const GModelSpectralNodes &model)
Assignment operator.
const std::string & name(void) const
Return parameter name.
Random number generator class definition.
#define G_INTENSITY_GET
Abstract spectral model base class.
virtual std::string type(void) const
Return model type.
int size(void) const
Return number of parameters.
void insert(const int &index, const GEnergy &energy, const double &intensity)
Insert node.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1253
virtual GModelSpectral & operator=(const GModelSpectral &model)
Assignment operator.
void extend(const GModelSpectralNodes &nodes)
Append nodes from node function.
void write(GXmlElement &xml) const
Set or update parameter attributes in XML element.
Definition: GModelPar.cpp:351
std::vector< double > m_mc_max
Upper boundary for MC.
std::vector< double > m_old_energies
Old energies.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
std::vector< GModelPar * > m_pars
Parameter pointers.
XML element node class.
Definition: GXmlElement.hpp:48
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:264
Spectral model registry class definition.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
Random number generator class.
Definition: GRan.hpp:44
virtual double flux(const GEnergy &emin, const GEnergy &emax) const
Returns model photon flux between emin, emax
virtual double eflux(const GEnergy &emin, const GEnergy &emax) const
Returns model energy flux between emin, emax
Spectral nodes model class definition.
std::vector< double > m_epivot
Power-law pivot energies.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
virtual int elements(void) const
Return number of GXMLElement children of node.
Definition: GXmlNode.cpp:586
Time class.
Definition: GTime.hpp:55
std::vector< double > m_mc_cum
Cumulative distribution.
Gammalib tools definition.
void reserve(const int &num)
Reserve space for nodes.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:587
GNodeArray m_lin_energies
Energy of nodes.
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
Energy container class.
Definition: GEnergies.hpp:60
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
Spectral nodes model class.
std::vector< double > m_log_values
log10(value) of nodes
std::vector< GModelPar > m_values
Node values.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
void set_flux_cache(void) const
Set flux computation cache.
#define G_WRITE
virtual std::string print(const GChatter &chatter=NORMAL) const
Print node function information.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
void append(const GEnergy &energy, const double &intensity)
Append node.
Model parameter class.
Definition: GModelPar.hpp:87
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:279
virtual GModelSpectralNodes * clone(void) const
Clone spectral nodes model.
void copy_members(const GModelSpectralNodes &model)
Copy class members.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
void free(void)
Free a parameter.
virtual void write(GXmlElement &xml) const
Write model into XML element.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const
Evaluate function.
Energy container class definition.
void fix(void)
Fix a parameter.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:170
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
GVector log(const GVector &vector)
Computes natural logarithm of vector elements.
Definition: GVector.cpp:1274
GModelSpectralNodes(void)
Void constructor.
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
virtual ~GModelSpectralNodes(void)
Destructor.
std::vector< double > m_mc_exp
Exponent for MC.
void update_pars(void)
Update parameter mapping.
void set_eval_cache(void) const
Set evaluation cache.
GChatter
Definition: GTypemaps.hpp:33
GEnergy m_mc_emax
Maximum energy.
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:235
GEnergy energy(const int &index) const
Return node energy.
std::vector< double > m_lin_values
Values of nodes.
#define G_INTENSITY_SET
const GModelSpectralNodes g_spectral_nodes_seed
Interface definition for the spectral model registry class.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:249
#define G_MC
#define G_ENERGY_GET
double intensity(const int &index) const
Return node intensity.
#define G_ERROR_GET
std::vector< double > m_mc_min
Lower boundary for MC.
void free_members(void)
Delete class members.
double error(const int &index) const
Return intensity error of node.
virtual GXmlElement * element(const int &index)
Return pointer to GXMLElement child.
Definition: GXmlNode.cpp:640
void update_eval_cache(void) const
Update evaluation cache.
virtual void clear(void)
Clear spectral nodes model.
void init_members(void)
Initialise class members.
std::vector< GModelPar > m_energies
Node energies.
void init_members(void)
Initialise class members.
const double MeV2erg
Definition: GTools.hpp:45
int nodes(void) const
Return number of nodes.
void read(const GXmlElement &xml)
Extract parameter attributes from XML element.
Definition: GModelPar.cpp:229
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:1422
std::vector< double > m_prefactor
Power-law normalisations.
GEnergy m_mc_emin
Minimum energy.
Exception handler interface definition.
virtual void read(const GXmlElement &xml)
Read model from XML element.
void mc_update(const GEnergy &emin, const GEnergy &emax) const
Set MC pre-computation cache.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
#define G_REMOVE
#define G_INSERT
virtual GXmlNode * append(const GXmlNode &node)
Append XML child node.
Definition: GXmlNode.cpp:287
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
#define G_READ
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
#define G_APPEND
void remove(const int &index)
Remove node.
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 set_cache(void) const
Set pre-computation cache.
#define G_ENERGY_SET
std::vector< double > m_old_values
Old values.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1295
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::vector< double > m_gamma
Power-law indices.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489
GNodeArray m_log_energies
log10(energy) of nodes
void xml_check_type(const std::string &origin, GXmlElement &xml, const std::string &type)
Checks the model type.
Definition: GTools.cpp:1819